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Abstract 



In this thesis, the Dyson-Schwinger - Bethe-Salpeter formalism is investi- 
gated and used to study the meson spectrum at zero temperature, as well 
as the chiral phase transition in finite-temperature QCD. 

First, the application of sophisticated matrix algorithms to the numer- 
ical solution of both the homogeneous Bethe-Salpeter equation (BSE) and 
the inhomogeneous vertex BSE is discussed, and the advantages of these 
methods are described in detail. 

Turning to the finite temperature formalism, the rainbow-truncated quark 
Dyson-Schwinger equation is used to investigate the impact of different forms 
of the effective interaction on the chiral transition temperature. A strong 
model dependence and no overall correlation of the value of the transition 
temperature to the strength of the interaction is found. Within one model, 
however, such a correlation exists and follows an expected pattern. 

In the context of the BSE at zero temperature, a representation of the 
inhomogeneous vertex BSE and the quark-antiquark propagator in terms of 
eigenvalues and eigenvectors of the homogeneous BSE is given. Using the 
rainbow-ladder truncation, this allows to establish a connection between 
the bound-state poles in the quark-antiquark propagator and the behavior 
of eigenvalues of the homogeneous BSE, leading to a new extrapolation tech- 
nique for meson masses. This is used to study the ground- and excited-state 
meson spectrum for all quark masses from light to bottom, for pseudoscalar, 
scalar, vector, axialvector and tensor mesons. Good agreement with exper- 
iment is found, e.g., for all ground states in the bottomonium system. 

In addition, new applications of the inhomogeneous vertex BSE, such 
as the possibility to calculate on-shell quantities like decay constants, are 
investigated. Finally, we study the influence of the infrared behavior of the 
effective interaction on properties of tt and p mesons. 
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This PhD-thesis has been submitted to the University of Graz in April 2011. 
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Chapter 1 

Introduction 



Within the standard model of elementary particle physics, the strong inter- 
action is described by Quantum Chromodynamics (QCD). This relativistic 
quantum field theory of quarks and gluons is a gauge theory with gauge 
group SU(3), such that the elementary degrees of freedom carry a quantum 
number termed 'color' [1], that can take three values. 

Some the most noteworthy features of QCD are asymptotic freedom [2, 3], 
dynamical breaking of chiral symmetry [4], and confinement. Asymptotic 
freedom means that at high energies the coupling of the theory becomes 
small, and perturbation theory is applicable. At low energies, where the cou- 
pling is large, QCD phenomenology is governed by the dynamical breaking 
of chiral symmetry which is responsible for the masses of the light hadrons. 
Hadrons, which appear in the framework of QCD as bound states of quarks, 
antiquarks, and gluons, are the only observed particles that interact strongly. 
This is usually referred to as confinement, and it is one of the most inter- 
esting problems in QCD, especially since an exact proof is still missing (for 
recent reviews, see [5, 6]). 

Due to confinement, the low-energy spectrum of QCD consists solely of 
bound states. This, together with the large coupling, necessitates the use 
of non-perturbative methods for any calculation of the physical spectrum 
of QCD. Two main branches have developed: Lattice QCD and functional 
approaches. 

In the case of Lattice QCD, a discretization of space-time leads to a 
regularization of the (Euclidean) path integral [7], such that observables 
like hadron masses can be calculated more or less directly. While these 
ab-initio calculations give a good description of the hadron spectrum (see, 
e.g., [8]), it is non-trivial to study the properties of the elementary degrees of 
freedom, quarks and gluons, in this approach. In particular, the theoretically 
interesting limits of very small and very large momenta are difficult to access. 

On the other hand, one may use a so-called functional approach and 
work directly with the Green functions of the elementary degrees of free- 
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dom. They appear as solutions to the Dyson-Schwinger [9, 10] or Functional 
Renormalization Group (FRG) [11-13] equations; if all Green functions were 
known, the theory would be solved. However, the Dyson-Schwinger equa- 
tions (DSEs) and FRG equations form an infinite tower of equations, and 
analytic statements about the solutions can only be made in special cases, 
e.g., in the limit of very small momenta (infrared analysis, see e.g. [14]). 

In this thesis, we use a functional approach to QCD based on Dyson- 
Schwinger equations, which are perfectly suited to study the properties of 
quarks and gluons in a nonperturbative framework. However, in order to 
provide a consistent description of QCD phenomenology, hadronic bound 
states have to be considered as well. They can be included via the so-called 
Bethe-Salpeter equation (BSE), which was proposed as a natural exten- 
sion of the Green function formalism to bound states [15-17]. As outlined 
in Chap. 2, the BSE and the DSEs can be derived consistently from the 
generating functionals of QCD, although it is interesting to note that the 
application of the BSE to (mesonic) bound states of the strong interaction 
[18] predates the formulation of QCD. 

Practical calculations, however, require sophisticated numerical meth- 
ods, and the solution strategies and algorithms used in subsequent chapters 
are discussed in Chap. 3. Chap. 4 is devoted to an application of the DSE 
formalism to the restoration of chiral symmetry at finite temperature, and 
in Chap. 5 we use the same framework to study meson phenomenology. 
Various aspects of the different versions (homogeneous and inhomogeneous) 
of the Bethe-Salpeter equation are discussed, which allow a numerical cal- 
culation of meson ground states as well as excitations. Furthermore, the 
influence of the infrared behavior of ghost- and gluon propagators on meson 
phenomenology is studied. Chap. 6 summarizes the results and provides an 
outlook. 

Part of the material contained in this thesis is available in the following 
articles: 

M. Blank and A. Krassnigg, The QCD chiral transition temperature in a 
Dyson- Schwinger- equation context, Phys. Rev. D 82, 034006 (2010) 

M. Blank, A. Krassnigg, and A. Maas, p meson, Bethe-Salpeter equation, 
and the far infrared, Phys. Rev. D 83, 034020 (2011) 

M. Blank and A. Krassnigg, Light-meson properties from the Bethe-Salpeter 
equation (2010), To appear in the proceedings of "Quark Confinement and 
the Hadron Spectrum IX", arXiv:1011.5772 

M. Blank and A. Krassnigg, Matrix algorithms for solving (in)homogeneous 
bound state equations, Comput. Phys. Commun. 182, 1391 (2011) 

A. Krassnigg and M. Blank, A covariant study of tensor mesons (2010), Ac- 
cepted for publication in Phys. Rev. D, arXiv:1011.6650 
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Chapter 2 

Formalism 



The aim of this thesis is to study quarks and mesons in the framework 
of nonperturbative QCD in a functional approach. The tools employed are 
the (quark) Dyson-Schwinger equation and the Bethe-Salpeter equations for 
bound states, which — although originally derived in different settings - 
can be treated in a consistent way, as demonstrated in the present chapter. 

In the following, we sketch how the Dyson-Schwinger equations (DSEs) 
for one-particle irreducible Green functions and the Bethe-Salpeter equation 
(BSE) for the propagator of a quark-antiquark system, via two-particle ir- 
reducibility, rigorously follow from the underlying quantum field theory of 
QCD. On the same footing, starting from the appropriate generating func- 
tionals, it is also possible to obtain a consistent truncation scheme of the 
quark Dyson-Schwinger equation and the BSE, thus allowing for numerical 
investigations while still preserving the most important symmetries of the 
theory. In this work, we distinguish three types of Bethe-Salpeter equations: 
the BSE for the quark-antiquark propagator; the homogeneous BSE; the in- 
homogeneous Bethe-Salpeter equation for a vertex that connects quark and 
antiquark to a color-singlet current (vertex BSE) . The derivation of the ho- 
mogeneous BSE and the vertex BSE from the BSE for the quark-antiquark 
propagator is demonstrated in Sees. 2.4 and 2.5, respectively. 

2.1 One- and two-particle irreducibility 

This section briefly reviews some of the basic points which are useful for 
the derivations given later on in this chapter. For a more complete and 
pedagogical account on generating functionals and the steps leading to one- 
particle irreducibility and also to Dyson-Schwinger equations, the reader 
may consult Refs. [19, 20]. For an introduction to two-particle irreducibility, 
see for example Ref. [21]. 

In order to introduce the concepts of one- and two-particle irreducibil- 
ity (often referred to as 1PI and 2PI) in the context of the path-integral 
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quantization of a field theory one has to start with its basic building blocks, 
the generating functionals. Formulated in Euclidean space, the generating 
functional of the full Green functions of QCD is 

Z[J] =N Jv<j> exp (SqcdW +4>-J) j (2-1) 

where <\> = {ip, tp, A^, c, c} represents all (fermion, gauge boson, and ghost) 
fields of the theory. We consider the theory in Euclidean space, and the 
action Sqcd[<P] of QCD is given by 

S QCD [ij,i;,A,c,c} = J d 4 x ($($ + ig4 + rn)il> + \F^Fi»,+£gf) • (2-2) 

In contrast to lattice calculations, the functional approach considered here 
is gauge- fixed, which is reflected in the gauge- fixing Lagrangian C g f. For 
Landau gauge, it reads 

Lgf = ^( + * + 9f j %A^ j . (2.3) 

The gauge field A^ = A^tf is an element of the algebra su(3), and the t l 
denote the generators of the gauge group 577(3). The superscripts k 
denote the color components, and repeated indices are summed over. The 
field strength tensor F^ u appearing in the action is given by 

= (d,Al - d v A^ - u.r'-Aj.A 1 ;:) f . (2.4) 

The generating functional or partition function Z[J] is a functional of 
the sources J = {J^, J^, Ja, Jc, Jc}, and the scalar product <p - J indicates 
a sum over all types of fields, as well as summation/integration over all 
discrete and continuous indices and variables. The full Green functions of 
QCD follow from Z[J] via functional differentiation, for example a full two 
point function is given by 

z m [J] = 8 lm. ( 2 .5) 

6 J id J k 

Note that this expression still depends on the sources J; Since here only 
vacuum Green functions are considered, the sources have to be set to zero 
at the end of the calculation. 

To study only the connected Green functions of the theory, the definition 
of a new generating functional, 

W[J] = -Log[Z[J]}, (2.6) 

can be used. To get from connected to one-particle irreducible Green func- 
tions, one further step has to be taken, namely a Legendre transform 

T[<f>] = W[J] + <f>-J. (2.7) 
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It defines another generating functional, the 1PI effective action r [(f)], which 
is, in contrast to Z[J] and W[J], no longer a functional of the sources J, 
but of the classical fields 4> given by 

A* 6W[J] (0 K\ 

* = ~5T ■ (2 - 8) 

In a diagrammatic representation, a 1PI diagram is constructed such 
that it does not become disconnected by cutting one internal line. The 
1PI Green functions, however, are not a subset of the connected Green 
functions. This can easily be seen by the relation between the propagator 
G$ = 5 2 W[J]/(5Ji5J k ) and the 1PI function rjjf = 5 2 T[(f)]/ (5(f)i5cp k ), 

G { 2r { 3 = <% • (2.9) 

Thus, the 1PI function is the inverse of the propagator This is 

important when constructing a diagram: 1PI functions (which are sometimes 
referred to as (proper) vertex functions) always have to be connected to each 
other via propagators, i.e. they do not have external legs. 

The sources as well as the classical fields are defined as functions of one 
space-time point x, such that they may be called 'local'. This is well suited 
to study the properties and interactions of the elementary fields of the theory 
under investigation. Bound states in the context of a quantum field theory, 
however, are represented not by elementary fields but by composite fields. 

In order to study them, we follow [22] and introduce a new type of (bi- 
local) sources Jn~(x,y), which depend on two space-time points x and y 
and are able to connect two fields, as indicated by the indices i, j. In the 
presence of these sources the generating functional Z[J] reads 

Z[J] = N j V(f)exp(-SQc D [<f>] + J d 4 xd*y<Pi(x)Jik(x,y)My?) , (2-10) 

where we consider the case of vanishing local sources Jj = for simplicity. 
The indices i, k label the field type, and repeated indices are summed over. 
As before, the generating functionals of connected and irreducible Green 
functions can be defined as 

W [2] [J] = -Log[Z[J]], and (2.11) 

T [2] [B] = W [2] [J] + j d A xd i yB i] {x,y)J ]l {y,x) , (2.12) 

where Bij{x, y) denotes the classical bilocal field corresponding to the source 
Jji(y,x). T[B] is referred to as two particle irreducible (2PI) effective ac- 
tion, and functional derivatives of Tp] [B] are the 2PI Green functions. The 
functional derivatives of W[ 2 ] [J] are in this case not connected, but bi-locally 
connected. 
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In strict analogy to Eq. (2.9) in the 1PI case, one can show that the 
inverse of the bi-locally connected two point-pair function, which represents 
the propagator of a particle pair, is a 2PI function (see for example [21, 23]). 



2.2 Dyson-Schwinger and Bethe-Salpeter equations 

The Dyson-Schwinger equations (DSEs) provide relations among the Green 
functions of a quantum field theory, and are sometimes referred to as quan- 
tum equations of motion. They can be obtained from a master equa- 
tion, which in turn is derived from the generating functional Z[J] (see 
e.g. [19, 20, 24]). The master equation in the 1PI case reads 



ST[(f>] SSqcd 



+ 



5 2 W[J] 5 
5J5J ~86 



. (2.13) 



In the second term, the expression in the square brackets is the argument 
of the functional derivative of the action, meaning that at any occurrence 

of a field 4>i the corresponding fa + JTJT^S^. ^ as ^° ^ e inserted (summa- 
tion over repeated indices is implied). Eq. (2.13) is the DSE for a one-point 
function, therefore to get higher n-point functions one has to apply subse- 
quent functional derivatives. Alternatively, the derivation can be done in a 
diagrammatic way and has been implemented as the Mathematica-package 
DoDSE, as described in [25]. As already noted in connection with the gen- 
erating functionals, for vacuum Green functions the classical fields have to 
be chosen such that the sources vanish, J = 0. 

In this work mainly one Dyson-Schwinger equation is used, namely the 
quark gap equation. Its pictorial representation is given in Fig. 2.1. In 
Euclidean momentum space, the equation reads [19] 

S(p)~ l = Z 2 (ip + Z m m) + g 2 Z 1F J ^- lfl S(q)F a u (p, q)D^{p - q) . (2.14) 

It is a nonlinear, inhomogeneous integral equation for the renormalized 
dressed quark propagator S(p) which can be represented by two dressing 
functions A(p 2 ) and B(p 2 ), such that 

Sip)" 1 = \ 1 t>A{p 2 ) + B{p 2 ) . (2.15) 



The first term on the left-hand side of Eq. (2.14) is the inverse bare prop- 
agator Sq 1 = \p + Z m m (multiplied by the renormalization constant Z%) 
including the bare quark mass m as well as the mass renormalization con- 
stant Z m . The second term is usually referred to as (proper) self-energy. It 
involves the bare quark-gluon vertex 4p7^ (A a denotes the Gell-Mann 517(3) 
color matrices), the full 1PI quark gluon vertex Z\p T"(p, q) (with the renor- 
malization constant Z±f), the gluon propagator D^ u (p — q) and again the 
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1 7^ _ 



— O— = — » + 

5(p) 5n(p) 5(9) 

Figure 2.1: The quark DSE (2.14) in pictorial representation. The blobs 
denote dressed propagators and vertices. 

dressed quark propagator. The index a indicates color and is summed over. 
In the derivation of Eq. (2.14), the 1PI quark-gluon vertex appears as func- 
tional derivative of the effective action, <5 3 r[<^]/(<5-0 5A^ Sip), and the gluon 
propagator as second functional derivative of W[J], 5 2 W[J]/(5Ja^ SJa„)- 

In order to obtain a similar equation for the two-particle propagator of, 
e. g., the quark-antiquark system, := <5 2 W^ 2 ] [J]/{SJ^SJ^), one has 
to separate 'interacting' and 'non-interacting' parts. Following [22] (cf. also 
[21, App. A]), this can be done by writing the 2PI effective action T[ 2 ] [B] as 

r [2] [B] = Tr [S^B] - TrLog [B] + T[B] , (2.16) 

where T[B] represents the 'interacting contributions' to the 2PI effective 
action. The vacuum value of the classical field B^x, where the corresponding 
source vanishes, i.e., where = 0, is given by [26] 

Btf, = S , (2.17) 

where S is the exact quark propagator. In this case, the inverse two-particle 
propagator is 



/ G (2)W_ 62T [2]i B ] 



[2] ; ■ SB^5B^ 



From this equation, one can define the interaction kernel 

K - (2 19) 
such that Eq. (2.18) can be written in the more familiar form 

(G^y 1 = S^S" 1 - K , (2.20) 

or equivalently 

G$ =SS + SSKG ( *I , (2.21) 

which is the well-known Bethe-Salpeter equation (BSE) for the propagator 
of the quark-antiquark system. 

It should be noted that this construction is in a strict sense only valid if 
the quark and the antiquark have the same flavor; it is however possible to 
extend the procedure to the general case. 
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2.3 Truncation and models 

In principle, the quark gap equation and the BSE given in Eqs. (2.14) and 
(2.21), respectively, are exact equations. However, in practice the gap equa- 
tion is, via the hierarchy of Dyson-Schwinger equations, coupled to infinitely 
many other equations, starting with the DSE for the 1PI quark-gluon vertex 
and for the gluon propagator. The BSE contains the so-far unknown inter- 
action kernel K. Therefore, one has to resort to truncations of the system 
to obtain numerical solutions of these equations. There exists, however, a 
connection between the quark self energy and the interaction kernel, which 
may be used as a guiding principle and is derived in the following (cf. [26]). 

From the form of the 2PI effective action given in Eq. (2.16), it follows 
that 

M^i-B-' + fiS. (2 . 22) 

In the vacuum, where all sources vanish, one obtains from Eq. (2.12) 

sr [2] [B] 



X R - 

such that 



, (2.23) 



S- 1 = + . (2.24) 

This is an exact equation for the inverse quark propagator, and by compari- 

<5f f Si 

son with Eq. (2.14) it is clear that ss gives the self energy. Therefore, the 
self energy and the interaction kernel are (up to a sign) the first and second 
functional derivative of T[S], which opens a possibility to find a consistent 
and numerically feasible truncation of both equations. In addition, if T is 
chosen in accordance with the symmetries of the theory , these are manifest 
in the solutions of Eqs. (2.14) and (2.21), as discussed in Ref. [26]. This 
is especially important for the application of QCD to light hadrons, where 
chiral symmetry is essential. 

A popular choice of T is depicted in Fig. 2.2. It corresponds to the 
rainbow-ladder (RL) truncation of the quark self energy and the interac- 
tion kernel K and correctly implements chiral symmetry independent of the 
functional form of the effective interaction Q. This truncation is used for the 
majority of calculations presented in this thesis, and the different types of 
model functions and the corresponding parameters used are given in App. B. 



2.4 Homogeneous BSE and normalization 

The inhomogeneous BSE for the two-particle propagator, Eq. (2.21), has 
poles at those energies where the system exhibits a bound state. The corre- 
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Figure 2.2: Diagrammatic representation of T which leads to the well-known 
rainbow-ladder truncation, with the effective interaction Q (for details see 
App. B). 

sponding on-shell condition is 

Pf = -M? , (2.25) 

where the index i labels the bound state, whose total momentum and mass 
are given by Pi and Mi, respectively. 

(2) 

The propagator G)J may now be expressed in terms of pole contributions 
and regular terms (denoted by R(p,q,P)) as [18, 27] 

= E * ( "-ff!pr Pi) + *<*«.*■) ■ (2.26) 

i 1 

x{p, Pi) = S a (p + )r^ (p, Pi)S b (p-) is usually referred to as the Bethe-Salpeter 
wave function consisting of the Bethe-Salpeter amplitude r^i (p,Pi) and the 
dressed (anti)quark propagators S a,b . Here, the indices a, b indicate the pos- 
sibly different flavors of the quark and the antiquark. The charge-conjugate 
wave function x is given by 

X (p, -Pi) = [C X (-P, -PilC- 1 ]* , (2.27) 

where the charge-conjugation matrix C = 7274 and the superscript t de- 
notes the matrix transpose. In general, uppercase Latin letters refer to total 
momenta (P and Pi), lowercase letters to relative momenta of quark and 
antiquark (q, p, and k), and lowercase letters with subscripts ± to the mo- 
menta of quarks and antiquarks themselves (p±, q±, and k±), which are 
defined as p± = p ± r)±P with the momentum partitioning parameters rj±, 
■q + + ?]_ = 1. For a complete review of the kinematics of the Bethe-Salpeter 
equation and the momenta involved the reader is referred to App. A, where 
all necessary definitions can be found. 

Inserting Eq. (2.26) into Eq. (2.21) and equating residues at the bound 
state pole P 2 = Pj yields 

T {h] {q,Pj) = J Kik^P^ik^T^P^ik-) • (2.28) 
k 
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r m (?,P) 




Figure 2.3: The homogeneous BSE, Eq. (2.28). 



This is the well-known homogeneous BSE for the on-shell Bethe-Salpeter 
amplitude (BSA) T^(p, Pj), which is valid only at Pj = —Mj. In pictorial 
form, it is given in Fig. 2.3. 

Eq. (2.28), however, only determines the homogeneous BSA up to a con- 
stant factor and an additional normalization condition is needed. The canon- 
ical norm follows from the two-particle propagator, by demanding that the 
residue at the bound state pole is equal to one, i.e., that the decomposition 
(2.26) holds exactly. There exist different ways to obtain the normalization 
condition, see for example [27] and references therein. Here, we consider one 
pole term in with an arbitrary residue rj (in the neighborhood of the 
pole, the other terms may be neglected), 



x(p,Pi)x(q,-Pj) 
p 2 -p? 



(2.29) 



or equivalently 



G$(p,q,P) 



^ + (Pi) 



x(p,Pi) x(q,-Pi) 

P, ~ {Pi)a 



(2.30) 



where the simple pole in the variable P^ has been made explicit. 

In general, if a function f(z) has a simple pole in z and it can be writ- 
ten as f(z) = g(z)/h(z), then the residue of / at a point zq is given by 
Res(f,zo) = g(zo)/h' (zq). This can be used to obtain the residue of G,J, 
since from Eq. (2.20) it follows that 



G 



(2) 

[2] 



[S^S- 1 - K] 



(2.31) 



Therefore, the residues of the pole in P^ on both sides of Eq. (2.29) are 



dP, 



S a (q + r 1+ P)S b (q-r ] _P) 5{p 



K(p,q,P) 



P=P t 



Ti 



2(iU 



x(p,Pi)x(q,-Pi) ■ (2.32) 
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Reshuffling the terms and taking the trace in a functional sense as well as 
in Dirac-space, one obtains 



2(Pi)/i 



Tr 



X(9,--Pi) 



d 

dPa 



S a (q + V+ P)S b (q 




q p 



(2.33) 



P=Pi 



If T{ = 1, this is the canonical normalization condition for the Bethe-Salpeter 
wave function \ [18, 27]. After expressing % in terms of T^j and propagators, 
and working out the derivatives (note that T^j is defined analogously to x), 
it can be converted to a condition for the BSA (see, e.g., [28]). As a result, 
the canonical norm can be defined as 



l 



TO* g 

2P 2 ap, 



Tr 



S" 1 ^ + i7+P)f -P i )5 ,6 (g-)r [/l] (g,P i ) 
- / S a {q + )r [h] (q,-P i )S b (q-ri-P)T [h] (q,Pd 



+ y y 5 a ( g+ )r w (g,-p J )5 ft (( / „)K(p,( Z ,p)5 a (p + )r [h] (p,p l )5 b (p_) 

9 V 



p=p, 

(2.34) 



where all dependence on P has been made explicit. Thus, to get a unit 
residue in Eq. (2.29), Ti = 1, Tr^i has to be replaced by T^/^fN. 



2.5 Inhomogeneous vertex BSE 

The homogeneous BSA represents an on-shell quark-antiquark-meson ver- 
tex. In order to obtain off-shell information without having to solve the full 
inhomogeneous BSE for the two-particle propagator, it is desirable to define 
a general vertex that connects quark and antiquark to a color singlet current. 

(2) 

This quantity inherits the pole structure of the two-particle propagator G,J 
and is restricted to one (mesonic) channel, i.e. one set of quantum numbers. 

We now derive a Bethe-Salpeter equation for this vertex, and start from 
the inhomogeneous BSE (2.21) by projecting one side onto a suitable (renor- 
malized) current Tq, such that 

eg r = sst + ssk eg r . (2.35) 
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K{k,q,P) 



Figure 2.4: The inhomogeneous vertex BSE, Eq. (2.39). 

To obtain the legless structure of a proper vertex analogous to the homoge- 
neous BSA, the external legs on the right are removed, 

S^S^G^To = T + K r . (2.36) 

We now define the inhomogeneous Bethe-Salpeter amplitude as 

r = S^S^G® r , (2.37) 

and end up with 

T = T + K S ST . (2.38) 

In less symbolic notation, with explicit momentum arguments and matrix 
structure in Dirac space taken into account properly, this equation reads 

r(p, P) = r (p, P) + J K(p, k, P)S a (k + )T(k, P)S b (k„) , (2.39) 
k 

which is the inhomogeneous vertex Bethe-Salpeter equation (see, e.g., [29, 
30]). In pictorial form, it is given in Fig. 2.4. 

Note that this definition of the inhomogeneous BSA is consistent with 
the form of the vertex BSE used in Ref. [29], which can be obtained by 
introducing the renormalized fully amputated quark-antiquark scattering 
matrix M as 

M = K + KSSK + ... = S~ 1 S- 1 G ( ^S- 1 S- 1 -S~ 1 S- 1 , (2.40) 
and inserting this definition into Eq. (2.36). The result is 

r = r + MSsr , (2.41) 

which is equivalent to [29, Eq. (19)]. 

2.6 Finite temperature 

So far, only the vacuum, i.e. zero temperature and density, has been con- 
sidered. The functional formalism used here, however, is also well suited 
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to study QCD at finite temperature. An extensive discussion of thermal 
quantum field theory is presented in Ref. [31]; Here, only the basic concepts 
and details of immediate importance for the practical calculations are given. 

One starts form the observation that the partition function as known 
from statistical mechanics, 

Z = r Tr[exp(-pH)] , (2.42) 

with the Hamiltonian H and the inverse temperature f$ = 1/T, can be 
written as a path integral. The time variable, however, becomes complex, 
and instead of 'real time' t one has to use 'Euclidean time' r = it, with 
re [0,0]. 

Thus, the formalism developed before is applicable, but time is compact. 
As a consequence, one has to introduce boundary conditions in time, such 
that bosonic fields become periodic and fermionic fields become antiperiodic. 
If the equations are transformed to momentum space, it is easy to see that 
the energies, which in Euclidean notation are the fourth component of the 
4-momentum, are discrete. These Matsubara frequencies are given by 

£l n = InixT for bosons, (2.43) 

u n = (2n + 1)ttT for fermions. (2.44) 

Therefore, the integral over the energy has to be replaced by a sum over the 
Matsubara frequencies. 

Due to the differences in the treatment of the fourth momentum compo- 
nent compared to the other three, the equations no longer appear in covari- 
ant form. Physically speaking, the system has been coupled to an external 
heat bath, and the direction of Euclidean time corresponds to the rest frame 
of the heat bath, defined by the four vector u». If this additional vector is 
taken into account, it is still possible to use a covariant notation to inves- 
tigate the implications of the finite temperature formalism on the structure 
of Lorentz and Dirac tensors appearing in the equations (cf. [32, 33]). 

In this thesis, the main object under investigation at finite tempera- 
ture is the quark propagator, whose inverse at zero temperature is given in 
Eq. (2.15). In that case, it contains the tensor structures f and 1, which 
are the only Lorentz scalars that can be constructed from the four-momenta 
7 M and p^. At finite temperature, with the additional four- vector Uu, the 
available constructions are i> and (p -u)^ for the 'vector' part of the propa- 
gator, as well as 1 and (jp-u)jpifi for the 'scalar' and a possible 'tensor' part, 
respectively. Note that the factors of p -u are chosen for later convenience. 
This leads to the general form [34] 

S(p, u)~ x = \i> A(p 2 ,u 2 ,p -u) + ip -u C(p 2 ,u 2 ,p -u)+ 

1 B(p 2 , u 2 ,p - u) + p - u jf> i/l D(p 2 , u 2 , p -u) (2.45) 
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with the four scalar dressing functions A, B, C, D. Choosing 

UfM = (0,0,0,1) and p M = (£>„) (2.46) 
gives the standard form of the quark propagator at finite temperature, 

S(p, 0J n ) = if A(p 2 , (jj n ) + i7 4 w n C(p 2 , U) n ) 

+ lB(f,u n ) +f lA uj n D{f,uj n ) . (2.47) 

The dressing function D, however, is power-law suppressed in the UV [34] 
and does not contribute in all cases investigated here. At zero tempera- 
ture but nonvanishing chemical potential it vanishes exactly since the cor- 
responding tensor structure has the wrong transformation properties under 
time reversal [32]. 

Similar to the quark propagator, the gluon propagator acquires addi- 
tional structures as well, such that it reads [31, 34] 

D^(k,n) = p^(k,n)G(k,n) + p^(k,n) F(k,n), (2.48) 

rp -> 0, a and/or v = 4 



P£(*,Si) = a„„ - - Pj„. (2.50) 



where G(k, fi) and F(k, f2) denote its two dressing functions. 

In order to calculate the quark propagator the corresponding DSE has 
to be converted to the finite temperature formalism as well. Therefore, the 
Matsubara sum and three-dimensional integral have to be inserted for the 
usual four-dimensional integration (J — > T'^^_ 00 JS). Together with the 
correct form of the quark propagator, Eq. (2.47), the quark DSE at finite 
temperature can be written as (see, e.g. [34]) 

l/Z^A(f,u n ) + l/Z! (l5(pV n )+i 74 w n C(pV n )) = 
Zf/Zf if + Z 2 /Z 1 (i 74 w n + m) + 

T fl J 9 2 D^(k,nh u ^S(q,ui)T^(p,u n -,q,u;i) , (2.51) 
l=-oo f 

where k = p — q and O = oo n — uii are the gluon three-momentum and 
Matsubara frequency, respectively, and Z\, Zf, Z 2 , Z^ are the necessary 
renormalization constants [34] . They are obtained from the renormalization 
condition that for a renormalization scale £, and p 2 + w 2 = £ 2 , 

S -1 (p,w ) = il -P + «74^o + m(Q, (2.52) 

which is the finite temperature analog of the renormalization condition at 
zero temperature, Eq. (3.3). 
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Chapter 3 

Numerical techniques 



In the present chapter, we give an account of the numerical methods em- 
ployed in this thesis to solve the gap equation at zero temperature and the 
homogeneous as well as the inhomogeneous vertex BSE, whose derivation 
was discussed in Chap. 2. 

Our basic point is that the homogeneous and vertex BSE in a fully 
numerical setup correspond to matrix equations, an eigenvalue equation in 
the homogeneous and a simple linear system in the inhomogeneous case. 
From this point of view, it is natural to apply well-known and efficient 
algorithms that are able to solve such problems for general complex matrices, 
as also advocated in [35]. Even though in this work we only solve the BSEs 
in rainbow-ladder truncation, the methods discussed here are applicable in 
a more general sense: in more involved truncations as well as in systems 
involving more than two constituents, like baryons which are described by 
a covariant three-body equation, cf. [36, 37]. 

As a prerequisite, we discuss the solution of the quark propagator for 
real and complex arguments including the renormalization procedure used 
here. In Sec. 3.3, following the arguments given in [35], we turn to a de- 
tailed description of the numerical representation of the BSE kernel matrix, 
which is the basic building block of the homogeneous as well as of the vertex 
BSE, and discuss the application of the implicitly restarted Arnoldi factor- 
ization to the homogeneous and of the Bi-Conjugate Gradients stabilized 
(BiCGstab) algorithm to the vertex BSE. 

To illustrate their efficiency, we provide in Sec. 3.4 a detailed descrip- 
tion of the application of our methods to the case of pseudoscalar quantum 
numbers and a comparison to more conventional approaches. In Sec. 3.5, 
the comparison is extended to other quantum numbers, and in Sec. 3.6 we 
directly compare the homogeneous and inhomogeneous vertex BSE in terms 
of numerical efficiency. 



15 



CHAPTER 3. NUMERICAL TECHNIQUES 



3.1 Quark propagator and renormalization 

As already stated before, our main tool to investigate the properties of 
quarks is the quark DSE, which reads 

S^y 1 = Z 2 (ip + Z m m) + g 2 Z 1F j ^S(q)T a u (p, q)D fll/ (p - q) . 

q 

This equation is exact, but we can not solve it directly since the gluon 
propagator and the quark gluon vertex needed as input are not known 
a priori. Therefore, we rely in most of our calculations on the so-called 
rainbow-ladder (RL) truncation, cf. App. B. In this truncation (and using 
Z2 Z rn m = Z4 m(ij, 2 )), the gap equation becomes 



S{p) = Z 2 \f + Z 4 m(/i ) + J —1^ S{q) —lv _ ^ 2 T^ v {p - q) , 

Q 

(3.1) 

with T^ v (k) = Sfj, v denoting the transverse projector and Q{k 2 ) the 

effective interaction, which absorbs in this setting the dressing function of 
the gluon propagator, any momentum dependence of the quark-gluon vertex, 
the coupling g 2 , as well as the quark-gluon- vertex renormalization constant 
Z\p. Following [28], we denote by m(// 2 ) the renormalized quark mass at 
scale //. 

In order to solve Eq. (3.1), i.e. to obtain the two dressing functions A(p 2 ) 
and B(p 2 ) defined in Eq. (2.15), it is projected on the two tensor structures 
T\ = \i> and T 2 = 1 by using that Tr[Tj • Tj] = if i 7^ j. In total, 



Mp 2 ) z 2 + 3p2 J ^ _ J 2 ) ^2^+ 



(I 



q 2 A 2 (q 2 ) + B 2 (q 2 ) 
x( P - g + 2 ^-^ a - g) - g) ) (3.2a) 

Bip ) = Z A m{n ) + 4 / — . (3.2b) 

J ip-q) 2 q 2 A 2 (q 2 ) + B 2 (q 2 ) 

The most straightforward way to solve these coupled, nonlinear, inhomoge- 
neous integral equations numerically is direct iteration. First, we choose a 
parametrization of the momenta (cf. App. A), then we apply the so-called 
Nystrom or quadrature method (cf. [38, Chap. 4]), which amounts to re- 
placing an integral by a sum over suitable quadrature weights and points 
and neglecting the error term in order to discretize the integration variables. 
Choosing the integration points to discretize the momentum dependence on 
the left-hand side as well allows to iterate Eqs. (3.2). 
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Depending on the form of the effective coupling, the occurring integrals 
need to be regularized. This is implemented in a translationally-invariant 
form via a Pauli-Villars cutoff applied to the gluon momentum [39]. Sub- 
sequent renormalization then gives the renormalized dressing functions and 
the values of the renormalization constants Z 2 and Z^ . 

Following [28], we use a subtractive renormalization scheme and fix Z 2 
and by demanding that at a renormalization scale \i in the perturbative 
domain the propagator becomes bare, 

Aifj 2 ) = 1 and B(^) = m(/u 2 ) . (3.3) 

Writing the inverse propagator as 

S'\p 2 ) = Z 2 \f + Z 4 m(/i 2 ) + \f{A\p 2 ) - 1) + (B'(p 2 ) - m(fi 2 )) , (3.4) 



with the regularized but unrenormalized dressing functions (defined via 

Eqs. (3.2) if Z 2 = Z 4 = 1) denoted by A'(p 2 ) and B'(p 2 ), one obtains 
for nonzero quark mass the renormalization conditions 

A(p 2 ) = i + A'(p 2 )-A'(fi 2 ) (3.5a) 

B(p 2 ) = m(fi 2 ) + B'{p 2 ) - B'(fi 2 ) . (3.5b) 

The renormalization constants are then given by 



Z 2 = 2-A'{{i 2 ) (3.6a) 

gV) 

m(fj, 2 ) 



Z 4 = 2-^§. (3.6b) 



The chiral limit solutions are obtained for m(fj?) = 0, and in this case 

B(p 2 ) = B'(p 2 ) . (3.7) 

For our numerical calculations we follow Refs. [28, 40] and choose [i = 
19 GeV. 



3.2 Quark propagator for complex momenta 

In order to use the quark propagator as an input to a BS equation in Eu- 
clidean space it is not sufficient to know the form of the dressing functions 
for real arguments. It follows from the kinematics of the BSE that their 
arguments k± in general are complex vectors. As explained in App. A, the 
dressing functions A(p 2 ) and B(p 2 ) constituting the propagator have to be 
calculated on a parabolic region in the complex p 2 -plane which is defined 
by the total momentum of the quark-antiquark system. From Eqs. (3.2), 
it is clear that the quark propagator for any complex momentum p can be 
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obtained by just inserting the corresponding value and performing the in- 
tegration over the (real) momentum q, once the solution on the real axis is 
known. 

This method is straightforward, but some problems are encountered. 
First, the (numerical) integral over the loop momentum q only converges 
at reasonable rates as long as p is close to the real axis. Second, to com- 
pute the integral the effective interaction Q has to be evaluated for complex 
arguments, although it is in most cases designed to (effectively) reproduce 
the behavior of the running coupling of QCD, u s (p 2 ), only on the real axis 
(cf. [28]). 

A different approach has been proposed in [41]. It is based on a shift 
in the integration variable of Eq. (3.1), such that the external momentum 
appears in the quark and not in the gluon momentum and therefore Q(p 2 ) 
is only needed for real arguments. Iterating Eqs. (3.2) in this setup requires 
the knowledge of A(p 2 ) and B(p 2 ) not only on the real axis, but inside the 
whole parabola. 

If, however, A(p 2 ) and B(p 2 ) are assumed to be holomorphic in the region 
of interest, and one bears in mind that a holomorphic function is completely 
specified on a region by its values on a closed, surrounding contour, one can 
iterate on such a contour (ideally, one chooses a parabola which is cut at 
some point in the UV, cf. Fig. A. 2), and use e.g. the Cauchy integral formula 
to compute the values inside. The values outside the parabola in the UV 
needed for the iteration can be obtained by the straightforward analytic 
continuation of a fit of the solution on the real axis. 

The only remaining problem is the numerical stability of the Cauchy 
integral formula if points close to the parabola need to be calculated. Here, 
we apply the method proposed in [42], which was first used in this context 
in [43]. It is based on the Cauchy theorem, and allows to extract the values 
of a holomorphic function f(z) for any point zq inside a closed contour by 
numerical integration, such that J(zq) is given by 

/(--^(E«)/(Ef^). (3.S, 

where the w[zj] denote the quadrature weights corresponding to the points 
Zj on the contour. The sums run over all contour points j = l,...,n. 
Similar formulae exist for all derivatives of f(zo). In this work, in order to 
calculate the canonical norm of the homogeneous BSA, Eq. (2.34), we need 
only the first derivative of the dressing functions A(p 2 ) and B(p 2 ), which 
can be calculated using 

/w=(t«-/wE^)/(E^). 

(3.9) 



18 



CHAPTER 3. NUMERICAL TECHNIQUES 



Therefore, a numerical representation of the quark propagator in the 
complex plane is given by the values of A(p 2 ) and B(p 2 ) on points on a closed 
parabola (including the points themselves and the corresponding weights) 
plus a fit to the dressing functions on the real axis in the UV outside the 
parabola. Note that the evaluation of the propagator given in this form is 
numerically very efficient, since no interpolation is needed. 

3.3 Bethe-Salpeter equations 

Once the quark propagator is known in the complex plane, one can proceed 
to the numerical solution of the Bethe-Salpeter equations (2.39) and (2.28), 
as discussed in detail in Ref. [35] . 

The first step towards a numerical representation of these equations is 
the analysis of the Lorentz and Dirac structure of the respective amplitudes. 
This structure is a result of the particular representation of the symmetry 
properties of the state under consideration under the Lorentz group, includ- 
ing the state's parity and spin. Therefore, the (in) homogeneous BSA T(q, P) 
is decomposed into Lorentz-covariant parts Tj(g, P) and Lorentz-invariant 
parts F % {P 2 , q 2 ,q ■ P), respectively, as 

N 

F(q, P) = ^2 Ti(q, P) F l (P 2 ,q 2 , q -P) , (3.10) 

i=l 

where the number of terms N as well as the tensor structure of the Tj and T 
depend on the quantum numbers of the state (for the explicit construction, 
see App. C). The Ti(q,P), which carry the Lorentz and Dirac structure, 
are usually referred to as covariants, whereas we call the Lorentz and Dirac 
scalar quantities F l (P 2 ,q 2 ,q -P) the components of the amplitude F(q,P). 
The Ti(q, P) represent a basis for the BSA and one is to some extent free to 
choose the details thereof. 

For numerical convenience, the basis elements are constructed such that 
they are orthonormal with respect to the (generalized) scalar product 

Tr (Tj ■ Tj) = 5ij , (3.11) 

where any occurring Lorentz indices are understood to be summed over. 
The details of the construction for mesons of any parity and spin discussed 
in this thesis are given in App. C. Using the decomposition (3.10), the 
BSEs (2.39) and (2.28) can be rewritten as coupled integral equations of 
the components depending on the scalar products of the momenta via the 
corresponding projections on the basis Tj. 

Note, however, that the structure of the covariants Tj does not deter- 
mine all quantum numbers of the quark-antiquark system. For constituents 
of equal masses, the homogeneous amplitude is an eigenstate of the opera- 
tion of charge conjugation defined in Eq. (2.27), thus defining the quantum 
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number of charge-conjugation parity (C-parity), C = ±1. This determines 
the symmetry properties of the components with respect to qP (cf. App. D), 
such that in this case a restriction of the symmetry of the components allows 
to select either C = +1 or C = -1. Here, however, we consider a general 
dependence on q -P, such that the methods described in this section are also 
applicable in the case of constituents of unequal masses, where C-parity is 
not well-defined. 

We proceed by considering the integrand in the homogeneous and the 
inhomogeneous vertex BSE, which is in both cases given by 

K(k, q, P) S a (q, P) T(q, P) S b (q, P) . (3.12) 

The amplitude T(q,P) is expanded in the chosen Dirac basis Tj{q,P) and 
the result is projected on Ti(k, P). Doing so, one obtains a matrix structure 
in the space of covariants, and Eq. (3.12) can be written as a matrix- vector 
multiplication in this space involving the BSE kernel matrix IC l j(k, q, P), 

lC)(k,q,P)F3(P 2 ,q 2 ,q.P) = 



Tr 



Ti(k,P) K(k,q,P) S'faPftfaP) S b (q,P) {P 2 ,q 2 ,q -P) , (3.13) 



where the sum over the repeated index j is implied. 

The index j of the components F^(P 2 ,q 2 ,q -P) can thus be viewed as 
a vector index, which has to be contracted with the corresponding index of 
the kernel matrix /C* (fc, q, P). 

The next step is to make the dependence on the continuous momentum 
variables P 2 , q 2 , and q P numerically accessible. As in the case of the gap 
equation (cf. Sec. 3.1), we apply the quadrature method to discretize the 
integrals and use the same set of points also for the momentum dependence 
on the left-hand side. The homogeneous and the inhomogeneous vertex 
Bethe-Salpeter equations can then be written as matrix equations in the 
covariants and the discretized momenta and read 

in the homogeneous case, and 

F i,V = F ijP + )cf Q F^ Q (3.15) 

in the inhomogeneous case. The indices i, j label the components, the 
multi-indices V, Q stand for all discretized momentum variables (summa- 
tion over repeated indices is implied). The BSE kernel matrix K, 1 - g is the 
same in both equations, and subsumes the interaction kernel, the dressed 



propagators of the constituents, the Dirac- and Lorentz structure, and the 
discretized integrations. It is applied to a vector Ffj? or F % " V , representing 
the homogeneous or inhomogeneous BSA. 
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Figure 3.1: The five largest eigenvalues of the homogeneous BSE plotted 
over \J — P 2 . If A = 1 (indicated the dotted line), the bound state mass M 
is given by M = y/ — P 2 . The ground-state (leftmost intersection) solution 
vector has positive C-parity (pion), the second has negative (exotic) and the 
third again has positive C-parity (excited pion). 



3.3.1 Homogeneous BSE 

With the results of the preceding section, the homogeneous BSE, given in 
Eq. (3.14) in index notation, can be written as 

F[h] = & ■ F [h] (3.16) 

using matrix-vector notation. As already mentioned in Sec. 2.4, this equa- 
tion is only valid at the on-shell points of the bound states in the respective 
channel, i.e. at certain values of the total momentum squared P 2 = —M 2 , 
where n = 0, 1, 2, . . . numbers the ground- and all excited states in the chan- 
nel. To find such a value of P 2 , one investigates the spectrum of K, as a 
function of P 2 , since Eq. (3.16) corresponds to an eigenvalue equation (with 
the dependence on P 2 made explicit) 

X(P 2 )F [h] (P 2 ) = /C(P 2 ) • F [h] (P 2 ) , (3.17) 

where the eigenvalue A(P 2 ) = 1. In other words, to numerically approach 
a solution of the equation, a part of the result has to be already known, 
namely the values M 2 , or — more precisely — the mass of the state one is 
looking for. The way out is a self-consistency argument, where the eigenvalue 
spectrum is plotted as a function of P 2 and those points with A n (P 2 ) = 1 are 
identified: the largest eigenvalue determines the ground state, the smaller 
ones in succession the excitations of the system (cf. Fig. 3.1). 

A great variety of algorithms is available to numerically tackle these kinds 
of problems, and the most commonly used is a simple iterative method. 
Similar to the other algorithms discussed in this section, it relies on the 
multiplication of the matrix K, on a vector and can successively be applied 
to find also excited states by projecting on states already obtained, see 
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e.g., Ref. [44]. This simple method, however, is not able to resolve pairs of 
complex conjugate eigenvalues, which may occur in the BSE, cf. [45]. In 
addition its convergence properties are not favorable, as demonstrated in 
Sec. 3.4.2. 

These difficulties are overcome by the use of more advanced algorithms. 
For this purpose, we utilize the implicitly restarted Arnoldi factorization [46] 
(implemented as the MPI-based ARPACK library), which is also frequently 
applied in lattice QCD studies, e.g. [47]. 

3.3.2 Inhomogeneous vertex BSE 

In the most compact notation, the inhomogeneous BSE (2.39) can be written 

as 

F(P 2 ) = F {P 2 ) + IC(P 2 )- F(P 2 ) (3.18) 

where the matrix K{P 2 ) is identical to the one in Eq. (3.16), and the vector 
Fq is given by the decomposition of To according to Eq. (3.10), Tq = £\ TiFq 
together with the discretization of a possible momentum dependence. 

Again, the simplest method to treat this problem is a direct iteration. 
Mathematically, this corresponds to the representation of the solution by 
a von Neumann series (cf. [38, Chap. 4]), which can be shown to converge 
as long as the norm of the operator K, is smaller than one, ||/C|| < 1. For 
matrices, this norm can be related to the largest eigenvalue, such that for 
P 2 > —Mq, the iteration converges. When P 2 approaches the ground state 
position — Mq from above, the number of iterations necessarily grows, and 
no convergence is obtained if P 2 < —Mq, as demonstrated in Sec. 3.4.2. 

However, a solution is possible for any P 2 if one rewrites Eq. (3.18) as 

F = (1 - ICy 1 • F , (3.19) 

i.e., F is given by the inhomogeneous term Fq multiplied by the matrix 
inverse of (1 — JC). F can then be computed by e.g. inverting the matrix 
exactly, which has been successfully used to resolve bound-state poles in the 
inhomogeneous BSA, as shown in [30]. On the downside, the direct inversion 
of a matrix is computationally expensive, and it is not straightforward to 
parallelize the procedure. 

A better approach is to view Eq. (3.19) as a linear system whose solution 
is to be found. Equations of this type are very common and several algo- 
rithms have been developed for their solution. In particular, if the matrix 
(1 — /C) is big, Eq. (3.19) is a typical application for the so-called Conju- 
gate Gradient (CG) algorithms. Many types of these iterative Krylov-space 
methods are available. In the case of the BSEs considered here, the matri- 
ces involved are neither hermitian nor symmetric, such that a good choice 
is the well-known Bi-Conjugate-Gradients stabilized (BiCGstab) algorithm 
[48], which is widely used for example in lattice QCD (cf. [49, Chap. 6.2], 
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P 2 [GeV 2 ] 

Figure 3.2: Component F\(P 2 , 0, 0) of the inhomogeneous pseudoscalar am- 
plitude calculated using BiCGstab vs. the square of the total momentum 
P 2 . The vertical lines mark the pole positions, corresponding to the pion 
ground- and first excited state (J PC = H ). 



where also the algorithm is described in detail) . This is the method we used 
to obtain the results presented in Sec. 5.4. 

3.4 Illustration: the pseudoscalar BSEs 

As an illustration, we discuss in detail the application of the algorithms pre- 
sented above to solve the homogeneous and inhomogeneous BSEs for pseu- 
doscalar quantum numbers and compare their efficiency in terms of the num- 
ber of matrix-vector multiplications needed to achieve a specified accuracy. 
We employ the rainbow-ladder truncation, i.e. the rainbow approximation 
in the quark DSE together with a ladder truncation of the corresponding 
quark-antiquark BSE. As effective interaction we use the model proposed by 
Maris and Tandy [40] with parameter co = 0.4 and light quarks (for details 
concerning the truncation and the effective interaction, see App. B). 

3.4.1 Kernel setup 

Choosing the rest frame of the quark-antiquark system, and applying the 
parametrization and discretization described in App. A, the kernel matrix 
Eq. (3.13) in our setup becomes (we use the orthonormal pseudoscalar co- 
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Figure 3.3: The number of matrix-vector multiplications needed for conver- 
gence of the simple iteration (left panel) and the Arnoldi factorization (right 
panel), plotted against the number of eigenvalues computed at a (typical) 
fixed value of P 2 = —Mq for random ("rand"), optimized ("opt"), and ideal 
initial conditions. 



variants constructed in App. C) 



= - 3^8 »fe>W j*y g ^_ g )2 2) ^(p - ?) 

-1 

xTr [T i {p,P)^S(q + )T j (q,P)S(q-) lu ] , (3.20) 

where w[qf], w[z m ] denote the quadrature weights (T^ represents the trans- 
verse projector) and the replacements p 2 — > p 2 , z p — >■ z s , q 2 —tqf, z — > z m 
have been made in all occurring momenta to implement the discretization. 
Therefore, the indices i;j label the components and r,s;l,m the momen- 
tum space points. For these illustrative calculations, we use N q = 32 and 
N z = 24, such that K, has the dimensions (32,24,4) x (32,24,4). 



3.4.2 Numerical efficiency 

To compare the efficiency of our algorithm of choice to solve the homo- 
geneous BSE, the implicitly restarted Arnoldi factorization, to the stan- 
dard method we compute the one to four largest eigenvalues of /C and 
compare the convergence in terms of the number of iterations needed to 
obtain an accuracy of the eigenvector of e = 10~ 8 , at a (typical) value of 
P 2 = —Mq = 0.0527 GeV 2 . For the simple iteration, our stopping criterion 
demands that the absolute change in any element of the eigenvector from 
one iteration step to the next does not exceed the desired accuracy e, while 
the vector is normalized to one with respect to the standard scalar product 
of C n , in order to be comparable to the ARPACK-library which works with 
the same scalar product. 
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Figure 3.4: The number of matrix-vector multiplications needed for conver- 
gence of the simple iteration (circles) and the Arnoldi factorization (squares), 
plotted against the number of eigenvalues computed at a (typical) fixed value 
of P 2 = —Mq, for ideal initial conditions. 



We compare three different initial conditions for both iterative proce- 
dures. First, we use the default settings provided by ARPACK, which are 
based on a random vector (cf. the ARPACK users guide [50]) opposed to a 
truly random initial vector F? n of complex numbers (where real and imag- 
inary part are uniformly distributed in the interval [0,1]) for the simple 
iteration. This choice is denoted by "rand". Second, we use optimized ini- 
tial conditions which anticipate the UV behavior of the amplitudes while 
not imposing further symmetries by choosing 



which we denote by "opt". Third, we consider "ideal" initial conditions 
by utilizing the eigenvectors obtained at a nearby value of P 2 , here P 2 = 
-0.053361 GeV 2 , as if n . In ARPACK, if more than one eigenvalue is sought, 
the initial vector is set to the sum of the previously computed eigenvectors. 

The results are given in Figs. 3.3 and 3.4, where the sensitivity of both al- 
gorithms to the initial conditions and the comparison of the efficiency in the 
ideal case are shown. From Fig. 3.3 it is clear that the Arnoldi factorization 
is less sensitive to a change in initial conditions than the simple iteration, 
and in general is more efficient. Even in the ideal case (cf. Fig. 3.4) for the 
first eigenvalue the advanced algorithm is 36% more efficient (7 iterations 
compared to 11), and becomes even more advantageous for an increasing 
number of eigenvalues. 

A further interesting observation from Fig. 3.3 is that the Arnoldi factor- 
ization for random initial conditions was more efficient for three eigenvalues 
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than when only two were requested. This is most likely due to a "cluster- 
ing" of eigenvalues number two and three for the algorithm, an effect which 
appears for eigenvalues close together and is also related to the eigenvectors. 
In this particular case, eigenvectors two and three have opposite C-parity 
or z-symmetry (cf. App. D), which may make them more easily distinguish- 
able for the algorithm and more easy to obtain as a result. The ARPACK 
library is very efficient at evaluating all eigenvalues in such a cluster, while 
convergence is slower if one asks for only one or a few of the eigenvalues in 
the cluster (see also [50]). 

To solve the inhomogeneous vertex BSE (2.39) we apply the direct iter- 
ation (summation of the von Neumann series) and the inversion using the 
BiCGstab algorithm, in the setup described above for pseudoscalar quantum 
numbers. In this case not only the structure of the amplitude, but also the 
structure of the inhomogeneous term Tq determine the quantum numbers of 
the solution. Following [29] , a possible choice for pseudoscalars is 

T = ^475 , (3.22) 

with Z4 the renormalization constant from Eq. (3.6). With this choice (pseu- 
doscalar, positive C-parity), as described in App. D no poles corresponding 
to negative C-parity appear in the solution, as can be seen from Fig. 3.2. 
The curve shown in this figure has been obtained with the BiCGstab algo- 
rithm, because as described in Sec. 3.3.2 the direct iteration fails to converge 
if P 2 < -M 2 . 

This is demonstrated in Fig. 3.5, where the number of matrix-vector mul- 
tiplications needed for convergence is plotted against P 2 for both methods. 
Here, we compare two choices of initial conditions: "ideal" starting values 
analogous to those discussed for the homogeneous equation, where the result 
of the calculation for the preceding value of P 2 (starting from P 2 = 1 GeV 2 
in this case) is taken as initial guess, and the standard choice F; n = 0. For 
both algorithms, the inversion was calculated to an accuracy of e = 10 -8 , 
such that 

||F -(1-/C)-F fin || <e, 

where F^ n denotes the numerical result of the inversion, and the norm is a 
maximum-norm, which corresponds to the maximal absolute value of any 
element of the vector. 

It is clear that, independent of the initial condition, the number of 
matrix-vector multiplications needed for the direct iteration diverges as P 2 
approaches — Mq (note that Fig. 3.5 uses a logarithmic scale on the verti- 
cal axis). The inversion with BiCGstab, however, converges for all P 2 with 
nearly the same speed, needing approximately 10 matrix- vector multiplica- 
tions. 

A further observation from Fig. 3.5 is that very close to the pole the 
standard initial conditions lead to faster convergence (by 1 iteration) than 
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Figure 3.5: The number of matrix- vector multiplications needed for con- 
vergence of the iterative solution of Eq. (2.39) (circles) and the BiCGstab 
algorithm (squares), plotted on a logarithmic scale against the square of 
the total momentum P 2 of the amplitude, for "ideal" initial conditions (full 
symbols) as well as for Ff = (open symbols). The vertical line indicates 
the position of the ground state P 2 = — Mg of the system. Note that the 
straightforward iteration does not converge for P 2 < —Mq. 



the ideal ones. This is most likely connected to the sign change in the 
components due to the pole, which is not taken into account in our definition 
of ideal initial conditions. In addition, these initial conditions also strongly 
depend on the step size. 

3.5 Other quantum numbers 

In order to check whether the results described above are limited to the 
special case of pseudoscalar quantum numbers, we discuss the application 
of the algorithms presented in this work to states with different parity and 
spin. 

In the homogeneous case, we consider scalar, vector, axialvector and ten- 
sor quantum numbers and use the same model and parameters as above in 
the pseudoscalar case. The bound states in these channels have been inves- 
tigated thouroughly in [51] and recently [52], where all details concerning 
the parameter dependence and a comparison to experiment are given. 

For the different quantum numbers, we compute the first eigenvalue 
(which on-shell corresponds to the ground state) for the four values of 
J PC = ++ , 1 , l" 1 , 2 ++ at a typical value of P 2 employing the Arnoldi 
factorization as well as the simple iteration, for both optimal and ideal initial 
conditions, as explained above. Note that the number of momentum-space 
points used for the integration is the same in all cases, while the number of 
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Table 3.1: Number of matrix- vector multiplications needed to achieve con- 
vergence of the Arnoldi factorization (Arn.) and the simple iteration (Iter.) 
when applied to the homogeneous BSE for the indicated quantum num- 
bers J PC at typical values of the total momentum squared P 2 , for "op- 
timized" and "ideal" initial conditions (with starting values computed at 
-(V--P 2 -0.01) 2 in each case . 

Initial conditions: optimal ideal 





P 2 [GeV 2 ] 


Arn. 


Iter. 


Arn. 


Iter. 


0++ 


-0.509796 


15 


40 


10 


15 


1~ 


-0.599695 


15 


74 


10 


31 


1+- 


-0.739772 


20 


103 


10 


15 


2++ 


-1.254400 


20 


63 


15 


46 



covariants differs. For pseudoscalars and scalars, the basis consists of four 
covariants, while for the other quantum numbers eight covariants have to 
be used (see, e.g., [52]), which increases the size of the BSE kernel matrix. 

The results, collected in Tab. 3.1, show that in each case the Arnoldi 
factorization was more efficient. Indeed, the advanced algorithm seems to be 
even more advantageous when applied to more complicated systems, which 
advocates its use also in studies of e.g. baryons. 

3.6 Comparison of homogeneous and vertex BSE 

As shown in Sec. 5.4, it is possible to calculate observables like meson masses 
and decay constants by using only the inhomogeneous vertex BSE instead of 
the homogeneous BSE. Therefore, it is interesting to compare the numerical 
efficiency of the solution methods for both equations directly, in order to see 
if the use of the vertex BSE provides a numerical advantage. 

To this end, we compare the matrix vector multiplications needed to 
solve the inhomogeneous BSE using the BiCGstab algorithm for a wide 
range of total momentum-squared to the solution of the eigenvalue problem 
for the homogeneous equation (using the Arnoldi factorization) in the same 
range for one and three eigenvalues. This calculation allows to obtain the 
masses of the pseudoscalar ground state and the first excitation for positive 
C-parity in our setup. 

The results of the comparison are given in Fig. 3.6. Even though ideal 
initial conditions were used in all cases, the solution of the inhomogeneous 
BSE took fewer multiplications than the eigenvalue calculations in the re- 
gion of the on-shell points of the ground state (corresponding to the first 
eigenvalue) and the first excitation (corresponding to the third eigenvalue), 
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Figure 3.6: The number of matrix-vector multiplications needed for conver- 
gence of the inhomogeneous BSE using BiCGstab (solid line) compared to 
the calculations of one eigenvalue (dash-dotted line) and three eigenvalues 
(dashed line) with the Arnoldi factorization. The vertical lines mark the 
ground and excited state in the h channel at Pq = — 0.0527GeV 2 and 
Pi = -1.3315GeV 2 . 



which are indicated by the vertical lines in the figure. This shows that in this 
case it is more efficient in terms of matrix-vector multiplications to use the 
inhomogeneous BSE to calculate the mass spectrum, especially the excited 
state. 

Demanding an even closer numerical correspondence of homogeneous and 
vertex BSE, we have also tried to extract the entire homogeneous ground 
state amplitude from the inhomogeneous BSA by fitting the residues in P 2 
at each point of our relative-momentum-squared grid in a rather naive ap- 
proach. While this appears possible in principle, the limiting factors in this 
case are the accuracies of the fit as well as the determination of the residues 
at the pole of the inhomogeneous vertex amplitude under investigation. In 
our pseudoscalar example, we were able to obtain the leading component 
of the homogeneous amplitude (which corresponds to the covariant T\ in 
Eq. (C.4)) with an accuracy of 2.6% when compared to the solution of the 
corresponding homogeneous BSE. The other components showed somewhat 
larger deviations such that this procedure cannot be performed without re- 
sorting to more involved fitting methods. Thus we conclude that, while 
obtaining the homogeneous BSA from the corresponding inhomogeneous 
vertex BSE is not impossible, with comparable numerical effort a determi- 
nation directly from the homogeneous BSE is much more accurate. 

The results for the 1 channel are presented in Fig. 3.7. For the ex- 
cited state, the advantages of the inhomogeneous equation are even more 
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Figure 3.7: The number of matrix-vector multiplications needed for conver- 
gence of the inhomogeneous BSE using BiCGstab (solid line) compared to 
the calculations of one eigenvalue (dash-dotted line) and three eigenvalues 
(dashed line) with the Arnoldi factorization. The vertical lines mark the 
ground and excited state in the 1 channel at Pq = — 0.5997GeV 2 and 
Pi = -1.0682GeV 2 . 



pronounced than for the pseudoscalar channel, whereas for the ground state 
the two equations are almost equivalent in terms of efficiency. 
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Chapter 4 

The chiral phase transition 



The framework of Dyson-Schwinger and Bethe-Salpeter equations provides 
a consistent approach to QCD phenomenology, which is very successful in 
describing not only the properties of the fundamental degrees of freedom, 
the quarks and gluons, but also mesons and baryons in the vacuum. In 
this setup, various aspects of chiral symmetry and its dynamical break- 
ing have been studied some time ago via the DSE of the quark propaga- 
tor (see, e.g. [53] and references therein). For the calculation of bound 
states, sophisticated studies in different truncations have been performed 
(see e.g. [36, 51, 54-56]). 

As discussed in Sec. 2.6, it is straightforward to extend the formalism 
to finite temperature and chemical potential, and several studies consider 
the restoration of chiral symmetry at finite temperature and the behavior of 
meson masses at and beyond the phase transition at zero chemical potential 
[57-63]. However, the situation is complicated by the additional variables 
introduced via the Matsubara formalism, such that most of these studies 
resort to simplifications that are already overcome in the zero temperature 
case. In addition, certain aspects of the phase diagram [64-69] including the 
deconfinement transition [70-72] have also been studied in this formalism. 

In the present chapter, which is based on Ref. [73], we focus on quarks at 
finite temperature and study the chiral phase transition via the properties 
of the quark propagator in rainbow truncation. It was shown some time 
ago that certain classes of models in the QCD gap equation in this trunca- 
tion yield a second-order phase transition with mean-field critical exponents 
[74, 75]. While the main emphasis of that investigation was the critical be- 
havior and universality class of the phase transition, it also provided the 
corresponding values for T c , which ranged between 120 and 174 MeV among 
the model interactions investigated. In somewhat different setups, using a 
separable interaction kernel in the gap equation [58, 60, 61, 76] or neglecting 
retardation effects [63], one obtains values in the range between 110 and 146 
MeV. 
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The pointwise behavior of all these interactions as functions of the gluon 
momentum is rather different, although in all cases it was determined by 
adjusting the relevant parameters to meson phenomenology at zero temper- 
ature. Therefore, it is interesting to ask whether there is any simple relation- 
ship between the momentum dependence of the interaction and the value 
of T c . Such a relationship could e.g. be analogous to that found recently in 
studies of meson properties using a particular form of model interaction with 
a one-parameter setup. There the value of the free parameter determines an 
effective range of the intermediate- and low-momentum parts of the interac- 
tion. While ground-state properties of pseudoscalar and vector mesons were 
unaffected by variations in the model parameter, masses of excitations of 
any kind (orbital or radial) showed a strong and systematic dependence on 
the model parameter [43, 51, 77-79] thus identifying excited-state properties 
as prime targets to study the particular details of the effective interaction 
in the DSE formalism, in particular in the nonperturbative regime (cf. also 
Sec. 5.2). In an analogous fashion, the goal here is to identify T c as a quantity 
with the same capability. 

4.1 Modeling at finite temperature 

The QCD gap equation at finite temperature, Eq. (2.51), retains the same 
form it has at zero temperature (cf. Fig. 2.1), such that for a solution the 
quark-gluon vertex as well as the gluon propagator are needed. Here, how- 
ever, we apply the rainbow-truncation, where, as stated in App. B, the gluon 
propagator and quark-gluon vertex are taken as bare and multiplied by an 
effective interaction in order to mimic the behavior of the full system. At 
finite temperature, the effective interaction is modified such that it follows 
the structure of the gluon propagator, Eq. (2.48), which has two dressing 
functions. Defining s := k 2 + Q 2 + m 2 with the Debeye mass m g [80], the 
effective dressing functions are given by 



In the following, we investigate the impact of different forms of the effec- 
tive interaction Q on the chiral phase transition. We consider the models 
proposed by Munczek and Nemirovsky (MN) [81], Maris and Roberts (MR) 
[28], Maris and Tandy (MT) [40], and Alkofer, Watson, and Weigel (AWW) 
[82]. The different functional forms of these effective interactions and the 
corresponding parameter sets are collected in App. B. 






s 




(4.2) 
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Figure 4.1: The order parameter Bq vs. temperature T, shown exemplarily 
for model MT2. 



4.2 Scaling analysis and critical exponents 

With all ingredients of the QCD gap equation at hand one can obtain a 
solution numerically and study the chiral phase transition temperature T c 
as well as the nature of the transition. The corresponding order parameter is 
the chiral condensate; Here, however, we use the equivalent order parameter 
Bq := B(0,luq), since it is easier to access and can be calculated more 
accurately. A typical behavior of the order parameter for a chiral-limit 
solution is shown in Fig. 4.1 and indicates a second order phase transition. 
Indeed, it has been found previously that the rainbow-truncated quark DSE 
generally yields a second order chiral phase transition with mean-field critical 
exponents [74, 75]. Using the scaling analysis described below, we have 
confirmed this behavior for all interactions considered here and thus enlarged 
the set of interaction types it had been shown for by MT and AWW. 

In the case of a second order phase transition, the order parameter 
Bo(t, h) obeys the scaling laws 

B (t,h) oc {-tf *->0~ and (4.3) 

h=0 

B (t,h) oc h 1 ' 5 h^0+. (4.4) 

Here, t = — 1 is the reduced temperature and h = y the reduced mass, 
a measure for the explicit breaking of chiral symmetry by a non- vanishing 
current quark mass m. (3 and 5 are the critical exponents of the phase 
transition. 

Although it is straightforward to use Eqs. (4.3)-(4.4) to obtain f3 and 
5, this procedure requires to solve the gap equation at h = or T = T c . 
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Figure 4.2: Chiral susceptibilities xt(T) (left) and Xh(T) (right) for the 
MT1 model plotted versus temperature T for various values of h. The 
curves between the data points were obtained by interpolation using cubic 
splines. 



This is numerically difficult, such that a direct evaluation does not allow 
fits to extract T c with the necessary precision and reliably observe scaling. 
Therefore we exploit further scaling relations and use chiral susceptibilities 
for our analysis. They are defined by 



T 9B (T, h) 
XT ■= -T c 



Xh 



dT 
dB (T, h) 



, and (4.5) 

h fixed 



dh 



(4.6) 

T fixed 



The maxima of these quantities for nonvanishing h are referred to as pseud- 
ocritical points, x V t an d Xh '> respectively The corresponding pseudocritical 
temperatures are denoted by T^ 6 and T^ c . They are obtained as the max- 
ima of the chiral susceptibilities with respect to temperature, as depicted in 
Fig. 4.2. 

The pseudocritical points Xr° anc ^ X V h a ^ so ODe y scaling laws. Their 
behavior for T ~ T c and h ~ is described by 

a hr 1+1 ", (4.7) 
X p h c oc hM 1 -^ . (4.8) 

Following [83] we use these scaling relations to obtain f3, 5, and T c for all 
models and parameter sets given in Tab. B.l, as illustrated in Fig. 4.3. 



4.3 Chiral transition temperature 

As already mentioned above, all interactions yield a second order phase 
transition with mean-field critical exponents, regardless of the strength of the 
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Figure 4.3: Peak heights of the chiral susceptibilities Xr an< ^ plotted 
versus the reduced mass h on log-log scale, for the MT1 model. The lines 
are linear fits through the calculated points which are used to obtain the 
critical exponents according to Eqs. (4.7) - (4.8). 

interaction, its range, or its pointwise behavior in any particular momentum 
range. In particular, the presence of the 5-function term does not make a 
difference in this respect. Another general observation is the fact that the 
function D(p 2 ,ujk) in Eq. (2.45) is identically zero for all interactions in 
rainbow truncation. 

The values obtained for the chiral transition temperature T c , however, 
are rather different among the various interactions. In the following we will 
argue how this feature can be exploited to discern various forms within a 
given truncation. The results for T c are summarized in Tab. 4.1 for the 
four interactions on a range of model parameters well-used in meson phe- 
nomenology. 

The MN interaction has to be treated somewhat separately, since the 
only free parameter in this case is D whose choice completely determines Q. 
In terms of plain numbers it is interesting to see that MN yields the highest 
value for T c , followed by MR. AWW and MT give the smallest numbers. 

For a given interaction, T c increases with oj, which is illustrated in panel 
(a) of Fig. 4.4. A straight-forward interpretation of this observation goes 
back to u> representing an inverse effective range of the interaction. In this 
picture, T c is higher for an interaction, for which this range is smaller, i.e. 
the major strength in the interaction comes from a region defined by the 
larger momentum scale u. However, this effect is obviously minor compared 
to the differences with regard to the form of the interaction, in particular 
the appearance of the <5-function term. 

To further illustrate the details of the differences between the various 
forms for the interaction as well as to highlight their characteristic features, 
we plot T c as a function of three more specific properties of each model, 
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Table 4.1: Results for the transition temperature T c rounded to MeV for 
the interactions and parameter sets defined in Tab. B.l. 



u [GeV] 


N/A 


0.3 


0.4 


0.5 


Model 


MN 








T c 


169 








Model 




AWW1 


AWW2 


AWW3 


T c 




82 


94 


101 


Model 




MT1 


MT2 


MT3 


T c 




82 


94 


96 


Model 




MR1 


MR2 


MR3 


T c 




120 


133 


144 



namely in panel (b), D in panel (c), and the chiral condensate in panel 
(d). While the attempt to quantitatively correlate T c globally (i.e. across 
the different interactions) to any of these fails, a clear structure is visible. 
Throughout all panels of Fig. 4.4 the AWW and MT results are rather close 
together, but clearly separate from MR, which again is clearly separate from 
MN. 

Both oj D and D investigated in the two middle panels have been inter- 
preted as an "integrated strength" of the interaction. In fact, the latter is 
proportional to the integral over d 4 q of the interaction at zero temperature 
for all interactions investigated here, if one leaves out the UV part J-jjv i n 
Eqs. (B.2) and (B.5). Fig. 4.4 in panel (c) clearly shows that there is no 
overall simple dependence of T c on D. In the MT interaction, a constant 
value for u D leads to (almost) unchanged masses and decay constants for 
ground-state pseudoscalar and vector mesons, for which reason lu D can be 
termed "integrated strength" in this case instead of D. Again, Fig. 4.4 in 
panel (b) indicates no simple dependence of T c on u: D. In both cases, the 
absence of a simple relation is exemplified already by the fact that three 
points with the same ui D, the characteristic feature for MT, in Fig. 4.4 (b) 
and three points with the same D, the characteristic feature for MR, in 
Fig. 4.4 (c) each correspond to three different values of T c . 

Finally, in Fig. 4.4 (d) we plot T c as a function of an order parameter 
of chiral symmetry breaking, the chiral condensate, evaluated at T = 0. 
In a very simple picture, one could argue that, given a certain form of the 
dependence of the order parameter on the temperature (see, e.g. Fig. 4.1), 
an enlargement of the value of the order parameter at T = would natu- 
rally lead to an increase in T c . However, panel (d) clearly shows that the 
interactions used in our setup are not that simple. 

Overall, these observations lead to the important possibility to establish 
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Figure 4.4: The different colors and symbols on the lines denote the corre- 
sponding interactions, namely: black (circles) - AWW; red (X) - MT; and 
blue (diamonds) - MR. The result for MN (single green square) is only shown 
in panels (c) and (d), since it does not have a value of uj associated with it. 
The insert in panel (d) enlarges the region of interest for the w-dependent 
interactions. 

distinct ranges of values for T c accessible for each type of interaction, which 
has the potential to rule out certain forms in a given truncation. We note 
here that various effects need to be taken into account beyond the simple 
setup presently used; we will attempt to quantify a certain part of the trun- 
cation effect below. At this point, however, we can still try to attribute the 
clear differences apparent from all panels in Fig. 4.4 to the structure of the 
interaction. The object of interest in this respect is the appearance and rel- 
ative strength of the 5-function term in Q, Eqs. (B.6) - (B.5): In the AWW 
and MT interactions, Eqs. (B.4) and (B.2), there is no such term. In the 
MR interaction, Eq. (B.5), it carries half the strength of the coupling in the 
sense of the integral J d 4 qQ ~ D described above. In the MN interaction, 
Eq. (B.6), it is the only term and obviously represents all of the interaction's 
strength. From this, one can say that a less pronounced <5-function term in 
the interaction will lead to lower values of T c . 

4.4 Conclusions 

We have investigated the impact of various forms of an effective quark-gluon 
interaction in the rainbow-truncated quark DSE in QCD at finite temper- 
ature on the critical temperature of the chiral phase transition. The es- 
tablished mean-field behavior of this transition in the present truncation is 
confirmed throughout. Regarding non- universal aspects of the transition, it 
is apparent that there is no simple overall relationship between the distribu- 
tion of the infrared strength of the interaction and the critical temperature, 
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although 1) a trend can be extracted identifying the part of the interaction 
strength carried by a 5- function term in some of the models as a qualitatively 
distinctive parameter, and 2) within the various models simple relations of 
T c and relevant model parameters do exist. 

4.5 Corrections to rainbow truncation 

While the rainbow truncation of the gap equation provides a simple setup 
for computations, one has to keep in mind that the necessary corrections 
in any given circumstance may change its results considerably. Therefore 
we discuss here some of the most important effects from which a change, in 
particular to the critical temperature, is to be expected. 

In [73] it was argued that from a phenomenological point of view one 
can try to quantify the contribution of the rainbow-truncated gap equation 
to the full (untruncated) result. Such an attempt has been made for meson 
properties at zero temperature in [84] and has been extended to the baryon 
sector and exemplified for several hadronic observables in [85] . The idea is to 
change the values of the available model parameters such that the resulting 
hadron properties, e.g. the rho-meson mass, are deliberately overestimated 
by the rainbow-ladder truncation result. In this way one expects corrections 
beyond this truncation to bring it to the experimental value. 

While we do not want to discuss the particular assumptions made in 
[84, 85], the effect described there is relevant for our present study, since we 
start with model parameters fixed to meson phenomenology at zero tem- 
perature. To get an estimate of the effects from a change of parameters 
at T = on T c , we simply adapt the two parameter sets referred to as 
"A" (which corresponds to our MT2 parameter set) and "B" (which has 
increased strength in the interaction to overestimate the rho-meson mass as 
discussed above) in Ref. [85] and add T c to the Tab. I given there. In this 
table, a systematic ratio of ~ 0.74 is found for value in A divided by value 
in B for several hadron properties (e.g. m p , f n , m^, mA, etc.). Our values 
for T c in this respect are 94 MeV for set A and 129 MeV for set B, which 
produces a ratio of 0.73 and thus fits perfectly into the picture. 

For the present discussion here the relevant point is that parameter 
changes inferred from the phenomenological estimation of corrections ex- 
pected beyond rainbow-ladder truncation are in agreement with investiga- 
tions at zero temperature. In particular, one obtains an increase of T c by 
about one third, depending on the details of the corrections assumed. 
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Mesons 



Strongly interacting particles, the so-called hadrons, can be divided into two 
classes, baryons and mesons, which can be distinguished since the number 
of baryons is a conserved quantum number. Thus, a meson is characterized 
as a hadron with baryon number zero. 

Mesons, in their simplest form, may be considered as bound states of a 
quark (baryon number ^) and an antiquark (baryon number —5), within the 
constituent quark model. These types of models, which are very successful in 
describing the hadron spectrum, provide a consistent classification scheme 
for mesonic states according to their flavor content via group theoretical 
arguments (see e.g., the review on the quark model in the Review of Particle 
Physics [86]). Nevertheless, the calculation of the meson spectrum in terms 
of constituent quarks has certain drawbacks. Especially the masses of the 
low-lying pseudoscalar meson octet (ir, K, n mesons) are not easily described 
by constituents with masses of approximately 300 MeV. 

The strong interaction, however, is well described by the relativistic 
quantum field theory of QCD, which was introduced already in Chap. 1. 
The elementary degrees of freedom of QCD, quarks and gluons, allow the 
classification of hadrons in the same way as the quark model. In addition, 
the properties of the light pseudoscalar mesons can be understood from a 
theoretical point of view: Their comparatively light masses are explained by 
their appearance as (pseudo) Goldstone bosons of the dynamically broken 
chiral symmetry of massless QCD. 

Thus, in order to provide a consistent description of mesons within the 
standard model of particle physics, it is necessary to treat them as com- 
posite states in the framework of QCD. Non-perturbative techniques are 
essential in this context, and e.g. lattice-regularized QCD is very successful 
in describing the hadron spectrum for ground states (see, e.g. [8, 87] and 
references therein) and excitations (see, e.g. [88, 89] and references therein). 
A different approach, as outlined in Chap. 2, is provided by the DSE/BSE 
formalism, which we employ in the following to study mesons as relativistic 
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bound states of QCD. 

In this context, mesonic states have already been investigated extensively 
using various approximations (see, e.g., [51, 55, 56] and references therein). 
For the most part the homogeneous BSE has been used, which provides good 
descriptions not only for the ground-state spectrum for various quantum 
numbers, but also for further mesonic properties such as decay constants, 
charge radii and electromagnetic form factors [79, 90, 91], as well as hadronic 
decays [92-94]. Some investigations employ the inhomogeneous vertex BSE 
as well, but mainly as a description for the quark-photon vertex in form 
factor calculations, or to calculate correlations at finite temperature [57]. 

However, the homogeneous BSE and the inhomogeneous vertex BSE, as 
discussed in Sec. 2.5, are intimately related. In Sec. 5.1, we further explore 
these relations and use them to gain insight into the structure of the Bethe- 
Salpeter equations. This allows to investigate in detail the consequences of 
the rainbow-ladder truncation, which leads (as shown in Sec. 5.2.3) to a new 
extrapolation technique to calculate the masses of bound states beyond the 
usually accessible region. It is applied in Sec. 5.3 to calculate the ground- 
and excited-state meson spectrum. 

The understanding of the structure of the BSEs also gives rise to new 
applications of the vertex BSE (Sec. 5.4). We show in Sec. 5.4.1 that bound- 
state properties such as decay constants can be calculated solely from the 
vertex BSE, and in Sec. 5.4.2 that it is possible to remove bound-state 
poles from the solution of the vertex BSE while keeping all non-resonant 
contributions. 

In addition, the DSE/BSE approach links the fundamental Green func- 
tions of QCD directly to hadronic observables. In the rainbow-ladder trun- 
cation, this link is provided by the effective interaction, which can be related 
to the running coupling and thus the ghost- and gluon propagator dressing 
functions. The infrared behavior of these dressing functions is currently de- 
bated, and in Sec. 5.5 we investigate the influence of the different available 
forms on the properties of bound states, especially the tt and p mesons. 
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5.1 Spectral representation of the BSEs 

In the numerical approach discussed in Sec. 3.3.1, the homogeneous BSE 
is treated as an eigenvalue equation via the representation of the BSE ker- 
nel and the propagators as matrices. This allows to obtain (numerical) 
eigenvalues and eigenvectors, which on-shell describe bound states and their 
properties. However, from a mathematical point of view each diagonalizable 
matrix is completely specified by its eigenvalues and eigenvectors. Therefore, 
even if they are not physical, off-shell eigenvalues and eigenvectors of the 
homogeneous BSE can be used to gain insights into the system under inves- 
tigation. In the present section, we explore the theoretical consequences of 
a representation of the homogeneous BSE via its eigensystem, which reveals 
further deep relations between the homogeneous BSE, the quark-antiquark 
propagator, and the vertex BSE. 



5.1.1 Notation 



For this study of the structure of the BSEs, we adopt an unusual notation. 
We denote a Bethe-Salpeter amplitude and its charge-conjugate by 



T(p,P) = \T(p,P)) and F(p, -P) = (T(p,P) 



(5.1) 



respectively, reminiscent of the bra-ket notation commonly used in quan- 
tum mechanics. A pair of dressed quark propagators is represented by the 
operator D(p,P), such that 

D(p,P) = S a ( P+ )®S b (p^ . (5.2) 

Applying D(p, P) on an amplitude gives a Bethe-Salpeter wave function, 

D(p,P) \T(p,P)) = S a (p + )T(p,P)S b (p.) = \ X (P,P)) ■ (5.3) 

Analogously, the charge-conjugate wave function is given by 

(X(P,P)\ = X(P,~P) = S b (p-)T(p,-P)S a (p + ) . (5.4) 

Note that the flavor indices a and b labeling the propagators are absorbed 
into the symbols used for the operator D(p, P), the wave function {x(p, P)\, 
and the amplitude \T(p,P)), in order not to overcomplicate the notation. 
An 'inner' product of a wave function with an amplitude corresponds to 



( X (p,P)\T(p,P)) =Tr 



J S b (p^)T(p, -P)S a ( P+ )T(p,P) 



(5.5) 



while the 'inner' product of two amplitudes or two wave functions is not 
well-defined. The homogeneous BSE (2.28) is constructed by applying the 
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quark-antiquark scattering kernel, represented by the operator K{p,q,P), 
to a wave function D{q,P) \T(q,P)), 

\T(p, P)) = K(p, q, P)D{q, P) \T{q, P)) . (5.6) 

Here, the argument q, which appears in the operator K as well as in the 
amplitude \T(q, P)), is integrated over. This is similar to the matrix notation 
used in Sec. 3.3, where the loop integral was represented as the multiplication 
of the BSE kernel matrix fC on a vector representing the amplitude. 

5.1.2 Homogeneous BSE 

Using the notation defined above, the homogeneous BSE (2.28) for a state 
(or excitation) i can be written as 

\Ti(p, Pi)) = K{p, q, Pi)D(q, Pi) \Vi{q, Pi)) . (5.7) 

As mentioned already in Sec. 3.3.1, this is an eigenvalue equation for the op- 
erator K D with eigenvalue Xi{Pf) = 1. The eigenvector is the homogeneous 
BSA. To find a numerical solution of Eq. (5.7), it is usually generalized to 
arbitrary P 2 , 

Xi(P 2 ) \Ti(p, P)) = K(p, q, P)D(q, P) \Ti(q, P)) , (5.8) 

where \Ti(p,P)) is a right-eigenvector of K D to the eigenvalue Xi(P 2 ). 

Although these eigenvectors can (in general) not be straight-forwardly 
interpreted with respect to physical states, they can be used to obtain insight 
into the structure of the Bethe-Salpeter equations. It was demonstrated in 
Sec. 3.3 that the linear operator K D (the BSE kernel matrix K, in the nu- 
merical setup) can be approximated to arbitrary precision by finite matrices, 
and thus has a discrete spectrum. Consequently, it can be expressed as a 
spectral sum via its eigenvalues and (left- and right-) eigenvectors. The 
right-eigenvectors, as stated above, are the solutions of Eq. (5.8). To find 
an equation for the left-eigenvectors, we consider the charge-conjugate of 
Eq. (5.8), 

(T t (p,P)\Xi(P 2 ) = (Ti(q,P)\D(q,P)K(q,p,P) . (5.9) 

On-shell, this equation gives the charge-conjugate BSA Ti(p, — Pi), but it 
is not an eigenvalue equation for K D. If, however, Eq. (5.9) is multiplied 
from the right by D(p, P), 

(T t (p,P)\D(p,P)X t (P 2 ) = {Ti{q,P)\D{q,P)K{q,p,P)D{p,P) , (5.10) 

we see that the left-eigenvector of K D to the eigenvalue Xi(P 2 ) is 

(Xi( P ,P)\ = (T l (p,P)\D(p,P) . (5.11) 
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With all these ingredients, the spectral representation of the operator 
K D can be expressed as a sum over the 'outer' product of all right- and 
left-eigenvectors multiplied by the corresponding eigenvalues, 

K(p, q, P)D(q, P) = £ A;(P 2 ) ^(p, P)) ( X i(q, P)\ ■ (5.12) 
i 

Here, the left- and right-eigenvectors are mutually orthogonal and normal- 
ized such that 

( Xj {q,P)\T i (q,P))=6 jti . (5.13) 

5.1.3 Quark-antiquark propagator and vertex BSE 

The operator K D is the basic building block of all Bethe-Salpeter equa- 
tions, and with the results of the preceding section, the quark-antiquark 
propagator as well as the inhomogeneous BSA can be expressed as a sum 
over eigenvalues and eigenvectors of the BSE kernel matrix. 

We first consider the BSE for the quark-antiquark propagator. From 
Eq. (2.20), it follows that 

G$(p,q,P) = D(p,P)(l-K(p,q,P)D(q,P))- 1 . (5.14) 

Inserting the spectral representation Eq. (5.12) into the above equation leads 
to 

G{| (p, q, P) = D(p, P) x _ w p2) P )> (Xi(q, P)\ 

i 

= E i _ A *(P2) \Xi(p,P)) (Xi(q,P)\ ■ (5-15) 
If all possible channels are considered, this is an exact representation of 

(2) 

GrL (p, q, P) which shows the intimate connection of the eigenvalues and 
eigenvectors of the BSE kernel matrix to the physical quark-antiquark prop- 
agator, when they are calculated off-shell. 

The inhomogeneous vertex BSE (2.39) can be written as 

\T{p, P)) = |r ) + K(p, q, P)D(q, P) \T(q, P)) , (5.16) 

with the inhomogeneous term represented by the ket-vector |To). Therefore, 
the inhomogeneous BSA is given by 

\T(p,P)) = (1 - K(p,q,P)D(q,P)y 1 \T ) , (5.17) 

as also discussed in Sec. 3.3.2. If the spectral representation Eq. (5.12) is 
inserted, we obtain 

\T(p,P)) = E 1 _ A 1 i(P 2 ) \ r *(P> P )) (Xi(q,P)\r ) , (5.18) 

which corresponds to a spectral representation of the inhomogeneous BSA. 
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5.1.4 Canonical Norm 

For convenience, the eigenvectors considered here are normalized according 
to the condition (5.13). If, however, the on-shell amplitudes are used to com- 
pute physical observables other than the mass, they have to be normalized 
canonically. The appropriate normalization condition, derived in Sec. 2.4, is 
based on the physical decomposition of the quark-antiquark propagator into 
bound state poles, which in the notation employed in this chapter reads 

0g» (P.?.-P) = £jjr !*(>■ A'! + ,, P) , (5.19) 

i 1 

where we denote the canonical norm of state i by Mi. In the vicinity of one 
pole i, the other pole- and regular terms become negligible, and inserting 
Eq. (5.15) we see that 

1 i / p\\ i ( p\i 1 \xi(p,Pi)) {xi(.q,Pi)\ /, on s 
Xi{p,P)) (xi{q,P)\ = — — ^ • ( 5 - 20 ) 



l-A.CP 2 ) 1 Mi P 2 -P { 

In summary, the canonical norm is given by 



M l= lim ( l , Xl{P l ] \ . (5.21) 

Using the rules of de l'Hospital, it can be expressed as 

M = p lim^ (--^HP 2 )) = ~K(Pf) , (5-22) 

such that the canonical norm of an amplitude normalized according to 
Eq. (5.13) is given by the slope of the eigenvalue curve at the on-shell point. 
With different notation, this form of the normalization condition has al- 
ready been derived in [95], and its validity independent of the truncation 
was demonstrated in [55]. 
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5.2 Extrapolation of meson masses 

In most studies to date the BSE approach was used to investigate the ground- 
state spectrum. While a study of excitations is possible in principle [43, 77, 
78, 96], two main problems appear: First, the analytic structure of the quark 
propagator limits the range accessible for a direct calculation, and either 
sophisticated methods for analytic continuation (e.g., [97]) or extrapolation 
techniques (e.g., [30]) need to be applied. Second, one can not exclude 
the appearance of so-called spurious or unphysical solutions (see, e.g., [45] 
for a collection of arguments), such that the identification of the physical 
excitations in a channel is unclear, especially for states of "exotic" C-parity. 

However, in order to attempt a solution to the second problem, the 
properties of the states in question need to be known. This in turn is only 
possible if the first problem has already been solved. In addition, a study of 
mesons of higher spins involves calculation of bound states of higher masses, 
such that even for ground states extrapolation techniques are needed. 

Based on the results obtained in Sec. 5.1, we therefore develop a new 
method for extrapolating the eigenvalues of the homogeneous BSE beyond 
the limits set by the complex conjugate poles in the quark propagator, which 
is subsequently applied to the excitations of the pseudoscalar, scalar, vector, 
axialvector and tensor channel. 

5.2.1 Consequences of the truncation 

So far, we have considered the general case, and no truncation has been 
applied. From this point of view, we can now investigate the consequences 
of assumptions used in a numerical setup. 

In the following, we consider a truncation constructed such that the 
quark-antiquark scattering kernel is independent of the total momentum, 



The rainbow-ladder truncation is a prominent example where the above 
condition holds, which entails that the only total-momentum dependence in 
the homogeneous BSE (5.7) lies in the two propagators, combined into the 
operator D(p,P) in the notation employed here. 

From Eq. (5.12), the interaction kernel can be written as a sum over 
(right-)eigenvectors of KD by amputating the propagators on the right- 
hand side, such that 



Since K is independent of P, we can express the interaction kernel by using 
eigenvectors and eigenvalues at any value of P, which will be used in the 
following to establish a connection between the poles in the quark-antiquark 
propagator and the eigenvalues Aj(P 2 ). 



K(p,q,P) = K(p,q) . 



(5.23) 




(5.24) 
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5.2.2 Behavior of eigenvalues 

As explained in Sec. 3.3.1, the on-shell points denning the masses of the 
bound states in the channel under consideration are obtained by investigat- 
ing the behavior of the eigenvalues with respect to the total momentum- 
squared, X(P 2 ). If A(P 2 ) = 1, the system exhibits a bound state. In addi- 
tion, from Eq. (5.22) we see that the canonical norm can be calculated from 
the derivative of this eigenvalue curve, again at the point where X(P 2 ) = 1. 
Therefore, it is interesting to investigate the consequences of Eq. (5.23) on 
the behavior of these curves. 

An eigenvalue curve which crosses one is directly related to a pole in the 

(2) (2) 

quark-antiquark propagator (p, q, P). Focusing on one pole i, G^' (p, q, P) 
can be expressed as the pole term plus regular terms denoted by Ri(p, q, P), 
which contain the non-resonant as well as the other pole contributions. We 
have 

= W + R,(p, q ,P) . (5.25) 



As above, we assume spectral normalization and denote the canonical norm 

,(2) 
r [2] 



(2) 

by Mi. On the other hand, via Eq. (5.15), G)J is given by 



ogfa«.^ = E '*^- ) l%& ,J ° )l . <«■*> 



l - HP 2 ) 

and it follows that 

\Xj(p,P))(x 3 (q,P)\ 



*i(p, ?> p )=£ ' - T-xm + kl ' (5 - 27) 



where 



5 _ \Xi(p,P)){Xi(q,P)\ 1 \Xi(p,Pi)){Xi(q,Pi)\ ( , 

Rl - i - UW) Jfi P^Tf ■ (5 ' 28) 

In order to link the physical pole in the quark-antiquark propagator to 
the eigenvalue curve, we consider 



Tv[K G] = Tr 



K(pt q) J_ ^MM^ + K(pf q)Ri{q ^ P) 



1 \Xi(q,Pd)(Xi( p,Pi)\ 

'Mi P^-Pf 

(5.29) 

where the trace is taken in the functional sense as well as in Dirac (and 
flavor) space, analogous to Eq. (5.5). In a diagrammatic language, this 
operation corresponds to joining the two ends of one diagram. 

In the first term of Eq. (5.29), we can make use of the generality of 
Eq. (5.24) with respect to the choice of P, which allows to express the 
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kernel via eigenvectors at P = Pi. In the second term, Eq. (5.24) is used 
directly, leading to 



Tr [KG] = Tr 



1 \Xi(q,Pi))(Xi(p,Pi)\ 



P 2 -P? 



+ ]T A,(P 2 ) llVfo P)> {Tj(q, P)| Ri(q,p, P) 



. (5.30) 



Evaluating the trace and using the orthogonality relation (5.13) gives 

^ [KG] = M % (P 2 -P?) + ? Aj(p2) (r ^' P)l \ T ^P)) ■ 

(5.31) 

Using Eq. (5.27), the second term in Eq. (5.31) becomes 

]T A,(P 2 ) (Tjiq, P)| R t (q,p, P) |r>, P)> = 

j 

£ A,(P 2 ) (T^g, P)| £ ^ P ; P) l { f^ P)l |r,(p, P)) 

+ ^A J (P 2 )(r i ( (7 ,p)|p J ( (7 ,p,p)|r j (p,P)) , (5.32) 



which by use of Eq. (5.13) reduces to 
However, if Eq. (5.15) is used directly to evaluate Tr[K G], we obtain 



A,(P 2 ) 



(P2) 



(5.34) 



Upon combining Eqs. (5.31), (5.33), and (5.34), we arrive at 

1302) = M Pf) +E A ^ p2 ) < r ; (9. p ) I ^ («. ^ p ) I r* (p, p)> • 

(5.35) 

The singular part of the operator Pj, defined in Eq. (5.28), vanishes by 
construction near P = Pj, such that the pole part dominates the behavior 
of Aj(P 2 ). If we assume this to hold also further away from the pole, the 
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MP 2 ) 




Figure 5.1: X(P 2 ) for quantum numbers h and light quarks. The lines 
joining the data points represent a linear fit, the gray bands are the corre- 
sponding statistical error bands. 



second term in Eq. (5.35) can be neglected and the behavior of the eigenvalue 
with respect to P 2 is determined via 

MP 2 ) 1 , s 

(5.36) 



1 - \(P 2 ) M . (pa _ p.2) ' 

If this holds in a realistic setup, it can be used to extrapolate masses of 
bound states beyond the region that is usually accessible, as discussed in 
the next section. 



5.2.3 Numerical results 

The methods advocated in Chap. 3 allow the reliable calculation of the 
biggest eigenvalues of the BSE kernel matrix, such that the behavior of 
Xi(P 2 ) for at least i = 1,...,5 can be investigated on numerically safe 
grounds. 

From Eq. (5.36) it follows that 

ll{p2) = x^) ~ 1 = Ml (p2 " ^ ' (5 - 37) 

which is linear in P 2 and goes through zero at P 2 = P 2 , i.e. on-shell, 
which makes it suitable for extrapolations, provided that the second term 
in Eq. (5.35) is indeed negligible. 

For our investigation we use the MT model, for parameter sets MT1, 
MT2, and MT3 (cf. App. B), for light, strange, charm, and bottom quarks 
with current quark masses as discussed in [51]. Starting with pseudoscalar 
quantum numbers, the behavior of X(P 2 ) is studied for eigenvalues one to 
three with positive C-parity for light quarks. The result is shown in Fig. 5.1, 
where the calculated values of X(P 2 ) are plotted together with a linear fit 
and the corresponding statistical error bands. As can be seen, the deviation 
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MP 2 ) 




Figure 5.2: A(P 2 ) for quantum numbers 1 and light quarks, analogous to 
Fig. 5.1. 

from the linear behavior is small, even away from the on-shell point which 
is given via A(-P 2 ) = 0. 

In Fig. 5.2 the same situation is given for quantum numbers 1 . In this 
case, differences from the linear approximation occur at small P 2 , however in 
the relevant area around A(P 2 ) = the linear fit represents the data points 
with good accuracy. In order to use all available data while still putting 
more emphasis on the region close to the on-shell point, each data point was 
weighted by 1/|A(P 2 )|. 

It should be noted that an inspection of all curves investigated in the 
following shows that Figs. 5.1 and 5.2 represent extreme cases, and the 
deviations from the linear behavior do in general not exceed the ones shown 
in Fig. 5.2. 

In order to test the reliability of this method for extrapolation, we use 
it to interpolate masses that can be obtained also from a direct calcula- 
tion, and apply the procedure to the quantum numbers J PC = h , ++ , 
1 ++ , l" 1 , 2 ++ , for quark masses from light to bottom. The result is given 
in Figs. 5.3 - 5.4, where the symbols denote the fit results, and the lines 
represent the solutions obtained from a direct calculation. The errorbars 
of the interpolation are in all cases consistent with the directly obtained 
values, which are in turn consistent with the results for the ground states 
given in [51, 52] and for the first excited states for pseudoscalar and vector 
quantum numbers in [43]. This indicates that the extrapolations performed 
are reliable, and that at least in the cases investigated here the systematics 
are under control, since deviations from the linear behavior result in larger 
errorbars of the fit. 
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Figure 5.3: —M 2 vs. w for the ground state, first and second excitation 
for J PC = h (panel (a)) and ++ (panel (b)). In panel (c), for J PC = 
1 , also the third excitation is shown. The black lines correspond to the 
results of the direct calculation where available, the symbols are the results 
of the interpolation/extrapolation. The errorbars correspond to the 95% 
confidence interval of the linear fit. 
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Figure 5.4: — M 2 for the ground state, first and second excitation for J PC = 
2 ++ (panel (a)), 1 ++ (panel (b)), and l" 1 (panel (c)), analogous to Fig. 5.3. 
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5.3 Ground- and excited-state spectroscopy 

The results of the previous section allow a stable extrapolation of meson 
masses beyond the limits given by the pole-pairs in the quark propagator. 
Thus, we are in the position to study the mass spectrum of quark-antiquark 
states in the isospin-symmetric limit. We consider the ground states and 
first two excitations for the quantum numbers J PC = ++ , 1 ++ , 1 H , 
2 ++ , and the ground state together with three excitation in the 1 channel, 
in the MT model (cf. App. B). 

5.3.1 Extrapolation in uj 

In order to compare our results to the experimental values, we utilize the 
dependence of our results on the model parameter ui while uj D is kept con- 
stant. For the ground states in the pseudoscalar, vector and scalar channel 
this dependence is rather weak [51]. On the other hand, the masses of the 
axialvector and tensor ground states change significantly with varying uj 
[51, 52]. This also happens in the case of radial excitations, as can be seen 
from the plots collected in App. E (cf. also [43]). In general, the overall 
picture indicates better agreement with experimental data for higher values 
of UJ. 

Although the MT model was originally designed to reproduce observables 
for light pseudoscalar and vector mesons, we investigate the possibility that 
this model with the accompanying rainbow-ladder truncation is actually well 
suited to describe heavy mesons (cf. [84, 98]). That picture is supported by 
recent studies in Coulomb gauge [99, 100], which show that the RL trunca- 
tion becomes exact in the limit of infinite quark masses. Thus, we choose 
the bottomonium system to fix the value for the model parameters. In 
particular, we focus on the quantum numbers 1 where high-precision ex- 
perimental data are available. Furthermore, as described in [51], the ground 
state of this channel was used to fix the bottom quark mass. 

Fig. 5.5 shows the ui dependence of the masses of the ground state and 
the first excitation (open symbols). Both show linear behavior, and we use a 
linear fit to extrapolate the mass of the first excitation to the experimental 
value. By this procedure, we fix the value of uj to 

uj = 0.572 GeV , (5.38) 

indicated by the star-symbol in the figure. 

It is important to note here that an extrapolation of the masses along the 
lines of Sec. 5.2.3 was only necessary for the excited state for uj = 0.5 GeV. A 
further note concerns the extrapolated value of the 1 ground state mass, 
which does not agree completely with experiment. This small discrepancy 
is acceptable at this stage, since the aim here is not to fine-tune the bottom 



52 



CHAPTER 5. MESONS 



M[GeV] 



10.0 



9.7 



9.8 



9.9 



9.6 




Y(2S) 



9.5 




Y(1S) 



9.4 



0.3 



0.4 



0.5 




Figure 5.5: w-dependence of the masses of the ground state and first excita- 
tion for the b-b J PC = 1 channel. The lines represent a linear fit through 
the calculated values for oj = 0.3 , 0.4 , 0.5 GeV. The first excitation is ex- 
trapolated to the experimental value at oj = 0.572 GeV, which is indicated 
by the star-symbol. 

quark mass, D, and oj fully but to provide an overall picture with reasonable 
effort. This is also justified a posteriori by the following discussion. 

In the vector channel, we consider a total of four states (ground state 
and three excitations), in the other channel we calculate three (ground state 
and two excitations). In order to obtain results at our fixed value of oj, we 
use the same procedure as described above and extrapolate linearly. The 
data and the corresponding linear fits are shown in Figs. E.l - E.6, collected 
in App. E. 

5.3.2 Results and discussion 

In Figs. 5.6 - 5.9 we present the final results of the extrapolations. In all 
cases, the errorbars correspond to the la confidence interval of the linear 
fit, and reflect the quality of the linear extrapolation combined with the 
uncertainties in the mass extrapolation. 

The spectrum of bb mesons, given in Fig. 5.6, shows good agreement 
with experiment for all ground states, which supports the conjecture that 
the rainbow-ladder truncation is well suited to describe heavy quarks in 
Landau gauge. In the vector channel the ground state was used to determine 
the quark mass and the first excitation to fix the value for oj, such that 
these states are fitted to the experimental values. In the other channels, 
the mass of the first excited state (where measured) is better described by 
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Figure 5.6: Results of the extrapolation for bottomonium (represented by 
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lines). The experimental value for the masses of hb(lP) and h^{2P) were 
taken from Refs. [101] and [102]. The experimental uncertainties are indi- 
cated with gray bands (where visible), and the errorbars of our results cor- 
respond to the la confidence interval of the fit in u, evaluated at u) = 0.572 
(cf. Fig. 5.5). 
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Figure 5.7: Our results for charmonium compared to experimental data [86], 
analogous to Fig. 5.6. The experimental value for the mass of h c (lP) was 
taken from Ref. [103]. Note that in the 1 channel the first and second 
excitation are degenerate within errorbars. 



54 



CHAPTER 5. MESONS 



Isovector 



M[GeV] 



2.5 - 



2.0 



1.5 - 



1.0 



0.5 



o.oh 



-7r(1800) ( 



tt(1300) 



-p(1700)° 
-p(1450)° 



-p(770) u 



•bi(2240) 



-ao(2020) 



- ai (1930) b!(1960) 

-*5-a,(1640)° <P 
-a (1450)° 1 ¥ 

3,(1260)° b!(1235)° 

-ao(980)° $ * 



-a 2 (1990) 
a 2 (1700)° 



a 2 (1320)° 



Figure 5.8: Results of the extrapolation in the isovector channel (light quark 
masses), compared to the experimental spectra [86], analogous to Fig. 5.6. 



the second excitation appearing in our calculation. However, it is unclear 
whether all excitations appearing as the solutions of a BSE correspond to 
physical excitations of the system under consideration. Hence, we have 
to allow for the possibility to discard certain so-called 'spurious' states in 
order to achieve agreement with experiment. In the vector channel, this 
possibly applies to the third and fourth excitation, which are degenerate 
within err or bars. 

In the charmonium system, Fig. 5.7, the pseudoscalar ground state and 
first excitation are in good agreement with the experimental data. As before, 
the vector ground state was fitted to the experimental value to fix the quark 
mass. The first and second excitation in our calculation are degenerate and 
are both consistent with the ip(2S) meson, and the ^(35) state is described 
by our third excitation. The situation is similar for J PC = 2 ++ . In the 
case of J PC = ++ , 1 ++ , and l" 1 , the ground states show deviations from 
the measured values. This may indicate that due to the lighter quark mass 
the rainbow-ladder truncation is less adequate to describe these states, and 
corrections have to be taken into account. 

This picture also holds for light isovector states, as can be seen from 
Fig. 5.8. Especially in the case of psudoscalar and vector mesons, the ground 
state and the first two excitations are in good agreement, if the first calcu- 
lated excitation in the vector channel is discarded. 

Finally, our results for the light isoscalar channel are shown in Fig. 5.9. 
Since we do not consider mixing and work in RL truncation in an isospin- 
symmetric setup, this channel contains pure nfi states (which are degenerate 
with the corresponding nn isovector states) and the pure ss bound states, 
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Figure 5.9: Results of the extrapolation in the isoscalar channel, compared 
to the experimental spectra [86] (black lines), analogous to Fig. 5.6. The 
blue symbols (o, V, o, >) correspond to ss bound states, the red symbols (□, 
A, <, < s?) to nil bound states in our calculation. The experimental error of 
/o(600) is shaded in lighter gray, since it overlaps with the error of /o(980). 
In the case of the /2(1430) state, the dotted line is omitted, since in [86] 
the mass is given as M ~ 1430 MeV, and no experimental error is reported. 
Note that we do not consider mixing. 
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where n denotes the light (up/down) quark. In the figure, they are rep- 
resented by the red symbols (□, A, <, and the blue symbols (o, V, o, 
>), respectively. For J = 1 one has ideal SU(S)-&&vot mixing, and 
thus our approximation is trustworthy in this case. This also allows a con- 
sistent fit of the strange quark mass via the 0(1020). In order to achieve 
better agreement with the experimental data in the pseudoscalar channel, 
one would have to consider mixing effects and explicitly include contribu- 
tions from anomalous terms. However, such investigations are beyond the 
scope of this thesis. 
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Figure 5.10: The inverse of the first component of the inhomogeneous vector 
amplitude l/i ? i(P 2 , 0, 0), as a function of the square of the total momentum 
P 2 . The zero-crossings give the masses of the ground state and the first 
excitation. 

5.4 Applications of the vertex BSE 
5.4.1 On-shell quantities 

Mesons as bound states are described by pole terms in the quark-antiquark 
propagator, which allows the construction of the homogeneous BSE to study 
their on-shell properties, as shown in Sec. 2.4. However, these pole contri- 
butions also translate into the vertex BSE, such that the inhomogeneous 
BSA, like the quark-antiquark propagator, has poles at the positions of the 
bound states (cf. [104]), as can be seen by, e.g., comparing the spectral 
representations of the quark-antiquark propagator, Eq. (5.15), and of the 
inhomogeneous BSA, Eq. (5.18). Since the vertex BSE can be solved in a 
numerically efficient way [35], it is interesting to discuss methods to obtain 
on-shell quantities like masses and decay constants from the inhomogeneous 
BSA. 

The most straightforward way to calculate bound state masses in this 
setup is to locate the positions of the poles in the inhomogeneous BSA with 
respect to the square of the total momentum, P 2 . Since the poles appear in 
all components of the amplitude, this can be done by investigating the zero 
crossings of any component via an interpolation or a fit [30, 104]. 

As an example, we consider the inhomogeneous BSA for vector quantum 
numbers and light quarks using the parameter set MT2 (cf. App. B). The 
correctly renormalized inhomogeneous term |To) in this case is [91] 

|r ) = Z 2lfl , (5.39) 

with the renormalization constant Zi from Eq. (3.6). As discussed in App. D, 
this choice of |To) has no overlap with states of positive C-parity, such that 
only non-exotic poles appear. 
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Figure 5.11: Diagrammatic representation of Eq. (5.41). 

We calculate (cf. [104]) the bound state masses by fitting the inverse 
of the first component Pi(P 2 ,0, 0) of the inhomogeneous BSA as shown in 
Fig. 5.10, and obtain for the ground and first excited states 

p 2 = -0.5991 GeV 2 F? = -1.0692 GeV 2 (5.40) 

which is in agreement with the corresponding results from the homogeneous 
BSE presented in Sec. 5.2.3. 

In the BSE formalism, mesonic decay constants are related to the cou- 
pling /- of a (canonically normalized) on-shell homogeneous BSA to a suit- 
ably renormalized current (T\ , which in diagrammatic form is given in 
Fig. 5.11. In the notation employed in this chapter, the equation reads 

fi = -^(t\ Xi (p,Pi)) ■ (5.41) 

For pseudoscalar and vector mesons, (r| = ^27^75 and (r| = ^2/37^, re- 
spectively [29, 40, 105]. 

As described in [104], the projection (5.41) can be extracted also from the 
inhomogeneous BSA without resorting to any information from the homo- 
geneous amplitude. Starting from Eq. (5.18), we consider the corresponding 
projection /- (P 2 ) of the inhomogeneous BSA on the current (f|, 

/ W( P 2 ) = ( f \ D (P,p)\r(p,p)) =E i_ A 1 . (P 2 ) ( f lxifo p i)) Mi,Pj)\ro) ■ 

(5-42) 

The residue of f^ h \P 2 ) at the pole i is given by 

Res(/f HP^Pi) = ±-(r\ X i(p,Pi)} (Xi(g,P)|r ) = ^/ r /r ■ (5.43) 

Thus, in order to calculate /- one in addition has to determine the projection 
of the homogeneous BSA on the inhomogeneous term, /p . Using the same 
arguments as before, this quantity can be obtained via 

Res(/ r ; ?i) (P 2 ),P) = -^(r | Xj (p,P))( Xl ( (? ,P)|ro) = -^(/y 2 . (5.44) 
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Figure 5.12: The projection (P 2 ) for the vector channel for light quarks, 
defined in Eq. (5.42). The line represents a pole fit to the data points, where 
the pole positions were taken from the fit of Fig. 5.10. 



In total, 

fi = ResC/f (P 2 ), P0/ V^Cfrf (P 2 ), Pi) ■ (5.45) 
In the vector meson channel, the decay constant of the i-ih state is given 

by 

U = Pf/Mi , (5.46) 

where Mi = \J —Pf denotes the mass of the respective excitation, and 

r = ^2/37^, as stated before. Therefore, the inhomogeneous term and 
the current only differ by a factor of 1/3, and the decay constant can be 
calculated using solely the residues of f^ h \P 2 ) via 

f v = ^^Res^V^P,) • (5-47) 

These residues are extracted from the pole fit shown in Fig. 5.12, where the 
masses (5.40) are used as input. 

For the parameters used here, this gives the decay constants of the p- 
meson and its first radial excitation corresponding to the masses in (5.40). 
The results, as presented in [104], are 

f p = 0.213 GeV f p , = 0.079 GeV , (5.48) 

which agrees perfectly with the values obtained by using the homogeneous 
BSE. 
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5.4.2 Projection on bound-state poles 

The pole structure of the inhomogeneous BSA is not only determined by 
the chosen structure of the amplitudes, but also by the choice of the inho- 
mogeneous term \Tq). Therefore, only poles whose homogeneous BSA has 
nonvanishing overlap with \Tq) can appear in the solution of the vertex BSE, 
which explains why (for momentum partitioning rj = 1/2) only poles with 
the same C-parity as |To) contribute. 

This argument, however, can be turned around to construct an inhomo- 
geneous term that removes the contribution of certain poles, provided that 
the corresponding homogeneous BSAs are known. Consider, for example, 
the inhomogeneous term 

|r / (p)) = |r )-|r i (p,p i )) , (5.49) 

where |Tj(p, Pj)) is the homogeneous BSA of state i. Inserting it into the 
spectral representation of the inhomogeneous BSA, Eq. (5.18), gives 

|r'( P ,p)> = £ 1 _i. {p2) l r ;(^ p )> (xifeP)l r'o(g)) . (5.50) 

At P = Pi, the term corresponding to the ith pole cancels, and the inhomo- 
geneous BSA reduces to 

|r>,p)) = J2 1 _ * lijfoPO) ( Xi (g,P)|r ) , (5.51) 

which by construction has no contribution from the state represented by 
\Ti(q,Pi)), i.e. it is regular at P = Pj. 

For a demonstration, we use the parameter set MT2 with light quarks 
(cf. App. B) in the pseudoscalar channel, where |To) = Z475 [29] with Z4 
from Eq. (3.6). In Fig. 5.13, we show the result of this procedure when 
applied to the ground state, for the first component Pi corresponding to the 
pseudoscalar covariant T\ given in Eq. (C.4). As can be seen, the ground 
state is indeed removed from the inhomogeneous BSA (solid line) as com- 
pared to the standard calculation where all poles contribute (dashed line). 

If this procedure is applied to all bound states which can be calculated 
directly, it allows for more stable extrapolations of pole positions based 
on fits of the inhomogeneous BSA. In general, if more than one pole is 
considered, these types of fits become unstable and thus complicate the 
study of excitations. In addition, the method proposed here allows a direct 
investigation of non-resonant contributions to the inhomogeneous BSA, since 
only the pole part is removed. 
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Figure 5.13: Component iq (P 2 , 0, 0) of the inhomogeneous pseudoscalar 
amplitude vs. the square of the total momentum P 2 , with (solid line) and 
without (dashed line) projection on the ground state. The gray vertical lines 
mark the pole positions of ground and first excited state. 

5.5 Homogeneous BSE and the far infrared 

The infrared behavior of the gluon and the ghost propagators in Landau 
gauge Yang-Mills theory is currently debated, and two types of solutions to 
these propagators which differ only in the far infrared have been proposed. 

It is expected that the form of the ghost and gluon propagator in the far 
infrared has very little impact on phenomenological quantities like meson 
properties. While this is consistent with the assumption that the different 
solutions are selected by a gauge choice [106], it is also supported by recent 
lattice calculations [107, 108], which find the gauge-invariant contributions 
to the p meson to be only about half a Fermi large. 

In the BSE setup, the effect of the far infrared on properties of 7r and 
p mesons has been investigated in [109], where the results and analysis de- 
scribed in the following have been presented. 

5.5.1 Effective Interaction 

For this study, the rainbow-ladder truncation of the quark DSE and the 
BSE is used. As discussed in Sec. 2.3 (cf. also App. B), this setup requires 
an effective interaction as an input, through which the Yang- Mills sector of 
the theory enters the equations. 

The form of the effective interaction can be based on the QCD running 
coupling a(k 2 ) [28, 40], where one important ingredient is the perturbative 
UV behavior. When viewed in this context, the rainbow-ladder truncation 
amounts to replacing the Lorentz- and Dirac-part of the product of the 
dressed quark gluon vertex and the dressed gluon propagator by their bare 
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Figure 5.14: The input running coupling (5.53), where SC denotes the scal- 
ing solution and DC the decoupling solutions. For more details, see [106]. 



counterparts multiplied by the (modeled) running coupling, such that 

If k k \ 
^ f i(p,q)Df lu (p-q)-^l^a(k 2 )-^l5 tlu , (5.52) 

where k = p — q, and r„(p, q) and D fjLU {p — q) are the full quark-gluon vertex 
and gluon propagator, respectively. 

In the context of Yang- Mills theory, the coupling can be defined as [14, 
106] 

a{k 2 ) = a{fi 2 )(G(k 2 ,fi 2 )) 2 Z(k 2 ,fi 2 ) (5.53) 

in the so-called miniMOM (minimal momentum subtraction) scheme [110]. 
Herein G and Z are the dimensionless ghost and gluon dressing functions, 
respectively, which are obtained from the scalar part of the propagators by 
multiplication with k 2 . For the scaling- type solution this coupling is infrared 
finite, while for the decoupling- type solutions it is infrared vanishing like k 2 
[106]. Note that this does not contradict the conjecture that both could be 
just gauge choices, since the running coupling defined in Eq. (5.53) is not 
gauge-invariant . 

The functional form of the input running coupling used here is taken from 
a fit [111] to the solutions given in [106]. It is shown in Fig. 5.14, where it can 
be seen that the coupling quantitatively depends on the choice of solution 
below a momentum scale of 1 GeV, and qualitatively below 100 MeV. In 
contrast to the scaling-type solution (denoted by SC), the decoupling-type 
solution is possibly not unique [106]. For this study, three representatives 
(denoted by DC 1 - DC 3) are selected, one of them (DC 3) being in rather 
good quantitative agreement with lattice results [106], in order to study also 
the dependence on different variants of decoupling-type behavior. 

On the other hand, truncation artifacts are introduced by using the 
bare quark-gluon vertex even if a(k 2 ) is given via solutions of the Yang- 
Mills system. In order to correct this mismatch and make contact with 
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phenomenology, we introduce the effective interaction 



G{k 2 ) = a{k 2 )F{k 2 ,uo,D) , 



(5.54) 



which replaces a(k 2 ) in Eq. (5.52). Here, 



F{k 2 ,u,D) = l + 4Tr 2 k 2 -,exp-^ 



(5.55) 



represents a non-trivial momentum dependence of the quark-gluon vertex. 
In total, this provides a two-parameter model similar to the Maris- Tandy 
model, Eq. (B.2) [40]. 

The definition (5.55) preserves the IR as well as the UV behavior of the 
running coupling a(k 2 ), such that the impact of the IR on meson prop- 
erties can be studied in a setup where chiral symmetry and its dynamical 
breaking are implemented correctly (cf. Sec. 2.3). Note that some additional 
results for the quark-gluon vertex are already available, which focus on its 
further tensor structures as well as its impact on bound-state calculations 
(see e.g. [112-116] and references therein). However, these are either specific 
to the scaling case [112], from the lattice for small volumes [113], or from 
more involved calculations using functional methods [114-116] and cannot 
be used in the same way as the gluon input for the comparison aimed at 
here. 

5.5.2 Quark propagator 

The quark propagator, consistently given as the solution of the rainbow- 
truncated gap equation (3.1), is a necessary input to the BSE. Its properties 
are encoded in two scalar dressing functions, which may be written as 



in terms of A(p 2 ) and B(p 2 ) defined in Eq. (2.15). Z g (p 2 ) and M(p 2 ) are 
referred to as wave-function dressing and mass function, respectively. 

In order to calculate these functions, the parameters of the effective 
interaction as well as the current quark mass m„, which acts as a driving term 
for the non-perturbative M in the quark DSE, need to be specified. Using 
the procedure described in Sec. 3.1, we have a value of m q = 0.0047 GeV at 
our renormalization point of (19 GeV) 2 and use two degenerate flavors. The 
parameters D and u are fixed by the pion properties as discussed below. 

The resulting quark dressing functions are shown in Fig. 5.15. They no 
longer show any qualitative difference when comparing results corresponding 
to the scaling- and decoupling- type inputs; in fact, only a small quantitative 
difference remains. The qualitative uniformity is a direct consequence of 
dynamical chiral symmetry breaking, since it makes the quark effectively 




(5.56) 
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Figure 5.15: The quark wave-function dressing (left) and mass function 
(right) for the four choices of parameters of the quark-gluon vertex fixed to 
pion properties, see Tab. 5.2. 
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Figure 5.16: The momentum dependence of the effective interaction (5.54) 
for the parameter values D and oj of Tab. 5.1 (left) and for the optimized 
values of Tab. 5.2 (right). 

decouple in the infrared and thus blind to the qualitative differences. The 
smallness of the quantitative differences among the solutions for Z q and M 
turns out to stem from the requirement of accurately reproducing the pion 
properties. 

5.5.3 Properties of n and p mesons 

We present results for masses and decay constants, where the latter are 
sensitive to the Bethe-Salpeter amplitude of the meson under consideration. 
The numbers are compiled in two tables to highlight two different aspects of 
the investigation. First, we keep the same vertex ansatz for all four different 
Yang-Mills input sets. We adjust the parameters of the vertex ansatz (5.55) 
such that the scaling solution reproduces the pion mass and decay constant 
adequately. The resulting values for D and u are 0.8 GeV 4 and 0.8 GeV, 
respectively. Then we compute the masses m and decay constants /, m n , /„-, 
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m n 


U 


m p 


f P 


Scaling 


0.1417 


0.0932 


0.7443 


0.1492 


Decoupling 1 


0.1420 


0.0956 


0.7641 


0.1533 


Decoupling 2 


0.1401 


0.0897 


0.7155 


0.1439 


Decoupling 3 


0.1389 


0.0805 


0.6321 


0.1239 


Experiment 


0.1396 


0.0922 


0.7755 


0.1527 



Table 5.1: Results for tt and p mass and decay constant after adjusting the 
parameters of the vertex ansatz (5.55) such that the scaling solution best 
reproduces the pion properties. All quantities are given in units of GeV; our 
results are accurate to less than half an MeV. 



m p , and f p with the same parameters also for the different decoupling- type 
cases; the corresponding results are listed in Tab. 5.1 in units of GeV. One 
can see that the more the input coupling differs from the scaling case, the 
larger the differences for f T , m p , and f p become, while the general picture 
is not spoiled at a qualitative level. The small variation in m n is due to 
the influence of the axial- vector Ward-Takahashi identity. This increase in 
discrepancy when only changing the gluonic input has also been observed in 
the literature [117, 118]. 

Next, we offer another possibility: irrespective of whether the different 
inputs correspond to different gauges or are just different truncation as- 
sumptions, it is to be expected that the quark-gluon vertex should change 
accordingly when the running coupling changes. To investigate this scenario, 
we have fitted the parameters D and oj for each input individually to repro- 
duce the pion parameters. The corresponding results are shown in Tab. 5.2. 
As one can clearly see, already a very moderate change in the parameters 
without excessive finetuning yields agreement of m n , f w , m p , and f p for all 
inputs on the percent level. Thus, any consistent change of the quark-gluon 
vertex and the gluonic sector permits an equally good description of m^, 
f n , m p , and f p . While not guaranteed, it is at least possible that a similar 
statement holds for arbitrary other observables as well; some confirmation 
of our observation along these lines has already been found in [55, 119]. 

It is interesting to check how the effective interaction (5.54) looks for 
the optimized fits. More precisely, we investigate the combined effect of 
the coupling and vertex dressing functions on different momentum scales. 
This is shown in Fig. 5.16. It is nicely visible that the four cases start to 
differ appreciably only below 1 GeV, but are qualitatively different in the far 
infrared. This observation supports the argument that the p has a gauge- 
invariant structure which exists at a level of 0.5 fm or below [107, 120], and 
is not too sensitive to larger scales. 
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D 




m n 


u 


m p 


f P 


Scaling 


0.800 


0.800 


0.1417 


0.0932 


0.7443 


0.1492 


Decoupling 1 


0.790 


0.810 


0.1414 


0.0932 


0.7443 


0.1496 


Decoupling 2 


0.820 


0.785 


0.1411 


0.0936 


0.7484 


0.1502 


Decoupling 3 


0.870 


0.760 


0.1407 


0.0935 


0.7485 


0.1514 



Table 5.2: As Tab. 5.1, but D and uj have now been adjusted for each input 
separately to fit the pion's mass and decay constant. All quantities are given 
in units of GeV except for D which has units of GeV 4 . 
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Chapter 6 

Summary and outlook 



The combined Dyson-Schwinger - Bethe-Salpeter formalism provides a con- 
sistent, nonperturbative, continuum approach to QCD phenomenology. This 
thesis focused on mainly two aspects: the chiral phase transition, which was 
investigated using the quark gap equation, and the meson spectrum studied 
by means of the (homogeneous) Bethe-Salpeter equation. 

For the numerical solution of the homogeneous and inhomogeneous ver- 
tex BSE, the application of sophisticated matrix algorithms was discussed. 
For the homogeneous BSE, which is generally solved as an eigenvalue equa- 
tion, a reliable method to compute several eigenvalues and eigenvectors was 
investigated, which facilitates studies of excited states in this approach. The 
inhomogeneous vertex BSE was treated as an inhomogeneous linear system, 
which could be solved via approximate matrix inversion. This method al- 
lowed to resolve not only the ground state pole in the inhomogeneous BSA, 
but also further poles corresponding to the various excitations. Both of 
these methods were shown to be numerically advantageous when compared 
to a more conventional approach. 

Inspired by the numerical methods, a spectral representation of the ho- 
mogeneous BSE, the inhomogeneous vertex BSE, and the BSE for the quark- 
antiquark propagator was found. Its theoretical consequences allowed to 
uncover the connections between these types of BSEs, and led to new appli- 
cations of the inhomogeneous vertex BSE, such as the possibility to calculate 
on-shell quantities like decay constants without resorting to the correspond- 
ing homogeneous BS amplitudes. 

In addition to these general considerations, more detailed studies of the 
properties of quarks and mesons were performed. For these, the rainbow- 
ladder truncation was applied, which provides a description consistent with 
the requirements of chiral symmetry. 

In the finite-temperature formalism, the chiral phase transition was stud- 
ied by means of the quark gap equation. In this context, several different 
forms of the effective interaction and their impact on the chiral phase transi- 
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tion were investigated. In all cases, the mean-field behavior of the transition 
in rainbow truncation was confirmed. The transition temperature shows a 
strong dependence on the particular model used, although no simple over- 
all relationship of its precise value to, e.g., the strength of the interaction 
in the infrared and mid-momentum regime could be established. However, 
within one model, a change in the parameters was found to affect the value 
of the chiral transition temperature in a similar way as hadron properties 
at vanishing temperature. 

In the investigation of meson properties, the rainbow-ladder truncation 
allowed to establish a clear connection between the bound-state poles ap- 
pearing in the quark-antiquark propagator and the behavior of eigenvalues of 
the homogeneous BSE. From this, a new extrapolation technique for meson 
masses could be motivated, which considerably enlarges the range of bound 
state masses accessible by the BSE approach. Combined with the numeri- 
cal methods proposed, the ground- and excited-state meson spectrum could 
be calculated for all quark masses from light to bottom, for pseudoscalar, 
scalar, vector, axialvector and tensor mesons. Especially for the bottomo- 
nium system, the results for all ground states showed good agreement with 
experiment, and also for lighter quark masses the pseudoscalar, vector and 
tensor ground states were well described in the rainbow-ladder approxima- 
tion used. A good description of scalar and axialvector states as well as 
a consistent identification of the physical excitations, however, remains a 
challenge. 

In the rainbow-ladder setup, the effective action can be used to investi- 
gate the impact of different treatments of the Yang- Mills system on meson 
phenomenology. It was shown that the behavior of the effective interaction 
in the infrared, which originates from different versions of the Landau- gauge 
ghost and gluon propagators, has little effect on the masses and decay con- 
stants of the 7T and p mesons. 

The results obtained in this thesis open up various possibilities for fu- 
ture research. The extrapolation technique, for example, allows a study 
of ground state mesons with higher spin, and therefore an investigation of 
Regge trajectories in the DSE/BSE approach. Another natural extension 
of this work is a detailed study of excited states, which might lead to a 
consistent identification of spurious excitations. In addition, the numeri- 
cal methods discussed here are not limited to either pure quark-antiquark 
systems or the rainbow-ladder truncation; they allow numerically efficient 
investigations of baryons, tetraquarks or even more complicated systems, as 
well as studies in more involved truncations. 
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Appendix A 

Euclidean momentum space, 
kinematics and integration 



A.l Euclidean metric 



As originally proposed by Wick [121], we consider the Bethe-Salpeter equa- 
tions and consequently also the quark DSE formulated in Euclidean momen- 
tum space. In this setup, the metric is given by = 5^ u and the Clifford 
algebra of the Dirac 7 matrices has to be modified accordingly. For the in 
this case hermitian 7 matrices we use the chiral representation which reads 
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/ -1 \ 
10 
10 

\ -1 / 

/ 1 \ 
0001 
1000 
\ 1 J 



7i 72 73 74 is diagonal. Using the slash-notation, 
we denote the scalar product of a four-momentum with the vector 7^ by 



(A.l) 



In this setup, the (timelike) total momentum of a quark-antiquark system 
P satisfies P 2 < 0, and the corresponding on-shell condition for bound states 
reads P 2 = — M 2 , where M > is the mass of the state. 
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Figure A.l: Momentum flow of the loop-diagram present in the homogeneous 
and vertex BSE (cf. [104]). 



A.2 Kinematics and parametrization 

Since we work in Euclidean space, it is convenient to use 4-dim. spherical 
coordinates, such that the momentum integration is written as 

oo 1 1 2vr 

-1 -10 

and a general four-momentum kn is parametrized as 

kfj, = Vk 2 (y/ 1 — z 2 \/l — y 2 sin(^), sj 1 — z 2 \/l — y 2 cos(0), y\/ 1 — z 2 , z^j . 

(A.3) 

The loop term in the homogeneous and vertex BSEs, Eqs. (2.39) and 
(2.28), is a simple triangle diagram, and we use the momentum flow as 
defined in Fig. A.l. Here, the total momentum P is given by the difference 
of the (anti-)quark momenta k± = k ± rj±P. The momentum partitioning 
parameters, denoted by rj±, satisfy r] + +T]^ = 1. The relative momentum k of 
the BSA in the loop is therefore given by k = r/_k + + r] + k^, and the relative 
momentum on the left hand side is denoted by q and defined analogously. 

Hence, the BSEs have three relevant momenta: P (total momentum), 
q (relative momentum), and k (loop momentum). Subsequently, using the 
4-dim. spherical coordinates introduced in Eq. (A.3), we choose P to be in 
the rest-frame of the system, 

P= (o,0,0,v A P2) . (A.4) 

The other momenta are chosen accordingly and read 

q= 7?(0,0,^l-z2,z g ) (A.5) 

and _ 

k = Vk? (o,x/l-z 2 \/l-y 2 ,yx/l-z 2 ,z^ . (A.6) 
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In this parametrization, the integration J d(j) is trivial. 

The components of the amplitude on the left hand side of the BSEs 
Eqs. (2.39) and (2.28) depend on the scalar products P-q = z q \fq 2 ~\/~P 2 , P 2 , 
and q 2 . Inside the integral, on the right hand side, the components depend 
on P - k = z^/k^V P 2 , P 2 , and k 2 . Thus, the BSE kernel matrix /C defined 
in Sec. 3.3 induces the following mapping on the momentum variables 

(q 2 ,z g )^(k 2 ,z) , (A.7) 

such that the integration J dy does not add a dimension to /C, although it 
is not trivial. 

The parametrization (A.4)-(A.6), together with the on-shell condition 
P 2 = —M 2 leads to an imaginary total momentum 

P=(0,0,0,iM) (A.8) 

in the homogeneous BSE and also in the vertex BSE as soon as P 2 < 0. As a 
consequence the (anti) quark momenta become complex, and the arguments 
of the dressing functions of the quark propagator A(k±) and B(k±) read 

k 2 z = k 2 - rf ± M 2 + 2iT]±MzVk 2 ~ . (A.9) 

This defines a parabolic region in the complex /cj.-plane on which the dress- 
ing functions are needed. The boundary of the region plotted in Fig. A. 2 is 
obtained by setting z = ±1. 

As described in Sec. 3.2, we solve the gap equation for complex momenta 
by iterating the equation on the boundary of the purple region in Fig. A. 2, 
while the blue region is covered by a fit of the solution on the real axis, 
which is analytically continued to complex values. The closed contour used 
in that procedure is given by the union of the curves 

c a (t) = -Ai + (Ai + A 2 ) t 2 - it 2 y /X 1 (X 1 + A 2 ) te[-l,l], (A.10) 

and 

c b (t) = A 2 + it2VAi(Ai + A 2 ) t€[-l,l]. (A.ll) 

Here, Ai = rj±M 2 gives the vertex of the parabola and A 2 defines the point in 
the UV where the parabola is cut. In our calculations, we use A 2 = 10 GeV 2 . 

The maximal size of the parabola (and thus the maximal bound state 
mass accessible) is in most cases (e.g. the MT model, cf. App. B) limited 
by the appearance of complex conjugate poles in the propagator dressing 
functions. Then, the contour method described in Sec. 3.2 is no longer 
applicable since it assumes the dressing functions to be holomorphic, and 
the integrals in the BSEs are divergent. 
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• (-J] ± 2 M\0) 
■ (0 , ±2j] ± 2 M 2 ) 











Re[k ± 2 ] 



Figure A. 2: Sketch of the region in the complex fc^_-plane spanned by the 
arguments of the propagator dressing functions. Purple: Bounded region 
used for solving the gap equation; Blue: region covered by a fit of the 
solutions on the real axis. 

A. 3 Quadrature rules 

After choosing a parametrization of the momentum integration, the next 
step is to discretize the momentum dependence. In this work, we straight- 
forwardly apply the quadrature method, and replace 




where w[qf], w[z m ], w[y n ] denote the quadrature weights and q 2 , z m , y n the 
corresponding nodes. The factors of k 2 /2 and ■s/l — z 2 are absorbed in the 
weights. 

However, in order to achieve good convergence of the integral J_ 1 dz, 
one has to choose a quadrature which takes the factor y/l — z 2 into account 
explicitly. For solving integrals like 

I = J dzy / l-z 2 f(z) , (A.13) 

we therefore make use of a generalization of the Clenshaw-Curtis quadrature 
[122] and expand the integrand f(z) into Chebyshev polynomials of the first 
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kind (denoted by T n (z)), such that 

oo 

/(*) = f+E a ^)- ( A - 14 ) 

n=l 

The T n (z) are given by 

T n (z) = cos(n arccos(z)) if — 1 < z < 1 (A. 15) 

for integer n > 0. Inserting into Eq. (A. 13) and defining z = cos(#) leads to 

^ oo J 

I=Y d6 sin 2 (#)T o (cos(0)) + J2a n / d# sin 2 (#)T n (cos(0)) 

n=1 

71" 7T 

OO 

<i# sin (#) + ^a n / d9 sin 2 {9) cos(n#) 



2 



n=l 





7T 7T 

= a --a 2 ^- (A.16) 

For a numerical approximation of the integral /, the series Eq. (A. 14) is 
truncated at finite order N, such that 

N 

f{z)*f( N \z)=Y J "a n T n {z), (A.17) 

n=0 

where indicates that the first and the last term in the sum have to be 
halved. In this case, f <yN \z) is a Chebyshev polynomial, and the coefficients 
a n are exactly given by [122] 

Upon inserting Eq. (A. 18) into Eq. (A.16), we arrive at 

f-£gV(-(£)) (!-(£)). (A-. 

The above formula gives a quadrature rule Yl,k=i w [ z k]f(zk) with 

z k = cos (jpj (A.20) 

and 

w[z fc ] = ( l-cos( ^X] . (A.21) 



\ N 
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This quadrature has two main advantages: first, its convergence properties 
are good since it takes the weight function \/l — z 2 into account explicitly; 
second, it is "nested", which means that it allows, without additional nu- 
merical effort, to estimate the errors by comparing the result for a given 
number of quadrature nodes with the corresponding result for half of the 
number of quadrature nodes. 
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Appendix B 

Models and parameters 



As mentioned in Sec. 2.3, the truncation of the DSE/BSE system that is 
used for most of our calculations is the rainbow-ladder truncation. The 
relevant ("interacting") part of the corresponding effective action is depicted 
in Fig. 2.2, and leads to the substitution 

IV(p, q)D IHt {p - q) -> 7„^P (<W " ^) , (B.l) 

where k = p — q, and r«(p, g) and D fJiU (p — q) are the full quark-gluon vertex 
and gluon propagator, respectively. Thus, in the quark DSE as well as the 
BSEs, the product of the full quark-gluon vertex and the gluon propagator 
is replaced by the corresponding bare quantities, which are then multiplied 
by the so-called effective interaction Q. 

In the present work we investigate and compare several different forms 
of this effective interaction. At zero temperature, for the studies of quark- 
antiquark systems presented in Chap. 5, we use the interaction defined by 
Maris and Tandy (referred to as MT, [40]) in addition to the scaling- and 
decoupling- type inputs discussed in Sec. 5.5 (cf. [106, 109]). In the MT- 
model, the effective interaction is given by 

^-D^e-W+Tvvtf), (B-2) 



where 



V (k 2 \ - 4tt 7m 7T F{k 2 ) 

Fuv{k } - l/21n[r + (l + P/A2 CD ) 2 ] " 



As given in [40], 

F(k 2 ) = [1 - exp(-k 2 /[Am 2 t ])}/k 2 . 

The additional parameters are mt = 0.5 GeV, r = e z - 1, N f = 4, A Q ^ D 
0.234 GeV, and 7m = 12/(33 - 2N f ). 
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Table B.l: Parameter sets D and co used for the different interactions MT 
[40], Eq. (B.2); AWW [82], Eq. (B.4); MR [28], Eq. (B.5); MN [81], Eq. (B.6); 
together with corresponding observables in T = meson studies, where 
available. Numbers are in GeV except for D (GeV 2 ) and the chiral con- 
densate, whose dimension is given explicitly. For the MT model we use the 
quark masses of Ref. [51]. 



Model 




D 


-(99)o 


m n 


m e 


MT1 


0.3 


1.24 


(.243 GeV) 3 


0.230 


0.780 


MT2 


0.4 


0.93 


(.242 GeV) 3 


0.230 


0.774 


MT3 


0.5 


0.744 


(.239 GeV) 3 


0.232 


0.762 


AWW1 


0.3 


1.47 


(.245 GeV) 3 


0.135 


0.745 


AWW2 


0.4 


1.152 


(.246 GeV) 3 


0.135 


0.748 


AWW3 


0.5 


1.0 


(.251 GeV) 3 


0.137 


0.758 


MR1 


0.3 


0.78 


(.241 GeV) 3 


0.139 




MR2 


0.4 


0.78 


(.250 GeV) 3 


0.139 




MR3 


0.5 


0.78 


(.255 GeV) 3 


0.139 




MN 




0.5618 


(.115 GeV) 3 


0.14 


0.77 



For our investigation of the chiral transition temperature (which is pre- 
sented in Ref. [73]), we compare the MT interaction to several Ansatze, 
which were proposed by Alkofer, Watson, and Weigel [82], Maris and Roberts 
[28], and Munczek and Nemirovsky [81], denoted by AWW, MR, and MN, 
respectively. At finite temperature, they are given by 



S(s) = 

AWW: = D%rse- sl ^ (B.4) 
MR: = D (^5 3 (k)5 k ^ + ^se- s ^ +T uv {s), (B.5) 

4-7T 3 

MN: = D —6 3 (k)6 k -i, , (B.6) 

where s := k 2 + Q 2 + m 2 , and the bosonic Matsubara frequency f2 of the 
gluons is given by the difference of the fermionic ones of the quarks, Vt = 
cok — uji, as mentioned in Sec. 4.1. The continuation of the MT model to 
finite temperature is straightforwardly implemented by replacing k 2 — > s. 
For the MR model, the function Fuv(k 2 ), given in Eq. (B.3), is the same as 
in the case of MT. 

At vanishing temperature, all of these Ansatze provide the correct amount 
of dynamical chiral symmetry breaking as well as quark confinement via the 
absence of a Lehmann representation for the dressed quark propagator. In 
addition, MT and MR produce the correct perturbative limit of the QCD 
running coupling, i.e. they preserve the one-loop renormalization-group be- 
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havior of QCD for solutions of the quark DSE. Note that only for MT and 
MR the integrals occurring in the gap equation need to be renormalized 
(cf. Sec. 3.1). 

All of the effective interactions introduced above contain the parameter 
D, which describes the strength of the interaction. The other parameter ui 
represents an effective inverse range of the interaction, which is not present 
in the MN model. For each model, the parameter sets used are given in 
Tab. B.l. 

In order to study mesons, the value of the current quark mass has to be 
fixed in addition to the parameters D and u. For the calculations presented 
in Chap. 5 using the MT model, we use the values of Ref. [51] which trans- 
late, at our renormalization scale of fi 2 = 19GeV 2 , to m u / d = O.OlGeV, 
m s = 0.0685GeV, m c = 0.828GeV, and m b = 3.71GeV for light (up/down), 
strange, charm, and bottom quarks, respectively. 
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Dirac structure of BS 
amplitudes 

The basic ingredients for the construction of the covariant structure of the 
Bethe-Salpeter amplitudes are the four-vectors occurring in the equation, 
namely 7„, q^, Pn- Out of those vectors, the four covariants that constitute 
the most general scalar amplitude read 

1, f , 0, i/2[f,f>]. (C.l) 

For computational convenience, we would like to construct covariants which 
satisfy the orthogonality relation 

Tr[Tj • Tj] = for i + j (C.2) 

where the product implies a summation over all Lorentz indices. To achieve 
this, we apply the Gram-Schmidt procedure to the covariants (C.l), and 
obtain 

Ti = l, T 2 = f , T 3 = $-f>^-, T A = i/2ty,f>]. (C.3) 

Note that these covariants describe a scalar quark-antiquark system, i.e., a 
state of positive parity. In order to obtain a basis for a pseudoscalar state, 
the above covariants have to be multiplied by 75, resulting in 

T= 75 _ 76(g) _ 75(^-^) r _ ji75(^-fo) 

1 2' 2 2 V / ^' 3 2 v /( y g 2 ' 4 V^V-ti 5 '?) 2 ' 

(C.4) 

where all Tj have been normalized to satisfy Eq. (C.2). The structure of the 
covariants, however, does not restrict the C-parity of the resulting state, as 
discussed in detail in App. D. 

To consider mesons of higher spin, one has to construct Lorentz-tensors 
of the appropriate rank which are symmetric, transverse in all open indices 
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and Lorentz-traceless [123, 124]. To add open indices in a way which auto- 
matically satisfies Eq. (C.2), we define 



T 

h 



Q,-P,^, (C5) 



such that the vectors {P^ , qj^ , 7j T } are orthogonal to each other. 
Constructing orthogonal (in the sense of Eq. (C.2)) covariants for a vector 
meson is now straightforward: The scalar covariants are multiplied by 
to give the first four, and by yjp- to give the second four covariants. By 
construction, they also satisfy the transversality condition, since qj^ as well 
as are orthogonal to P^. Note again that a basis of the axialvector 
states is given via a multiplication by 75. 

To proceed to the tensor meson (quantum numbers J p = 2 + ), we first 
need all symmetric combinations of q^ and r y^ T of rank two, namely [52] 

Mp V = ll T ql + qlll T and (C.7) 

= (C8) 

These quantities are by construction transverse in all indices as well as or- 
thogonal. To implement the traceless-condition, they have to be orthog- 

onalized with respect to the transversely projected metric tensor, gz v = 
p p 

Sfiv pT-- Finally, we arrive at 



^-tiu^k, (cio) 



which - by multiplication with the four scalar covariants - give a full orthog- 
onal set of covariants for a meson of spin 2 and positive parity. 

Note that in the numerical implementation all covariants are normalized 
according to % = Tij ^/Tr[Ti ■ Tj]. 
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If the quark and antiquark that constitute a meson have the same flavor, 
then the state has the quantum number of charge-conjugation parity (C- 
parity) in addition to spin and parity. In mathematical terms, the BSAs 
T(p,Pi) are eigenstates of the operation of charge-conjugation, defined via 

f(q, Pi) = [C T(-q, Pi) C" 1 ]' = Vc T(q, P) , (D.l) 

where C = 7274 is the charge-conjugation matrix and the superscript t 
denotes the matrix transpose. The eigenvalue r] c = ±1 gives the C-parity of 
the state. 

As mentioned in Sec. 3.3, the BSA can be decomposed into covariants 
and components (Eq. (3.10)), 

N 

T(q, Pi) = £ T n (q, P t ) F n {Pf, q\ q -P t ) . 

71=1 

If the covariants are constructed such that they have a definite C-parity £ n , 

= [C T n (-q, Pi) C" 1 ]' = ( n T n (q, P) , (D.2) 

the C-parity is manifest in the symmetry of the components with respect 
to q -Pi. To see this, we first insert the decomposition (3.10) and Eq. (D.2) 
into Eq. (D.l), 

N 

f(q, Pi) = ^ F n (Pi, q\ -q -P)[C T n (-q, P) C' 1 ] 1 

71=1 

N 

= F n {Pf,q\ -q -PKnTniq, P t ) . (D.3) 

71=1 

For the (anti) symmetry of the components with respect to qPi we introduce 
the notation 

F n {Pf, q\ -q -Pi) = ( n F n (Pl q\ q -Pi) , (D.4) 
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with = ±1, and arrive at 

N N 
r 1c Y,Tn{q,Pi)F n {P?,q\q-P l ) = Y,(Un)T n {q,P l )F n {P? 1 q 2 1 q-P i ) . 

n=l n=l 

(D.5) 

Thus, the C-parity r] c is given by CnCm which is the same for all n if the state 
under consideration is indeed an eigenstate of charge conjugation. Therefore, 
in order to produce a state with r] c = +1, the product Q n C, n has to be 
always positive, and if one covariant has negative C-parity, £ n = —1, the 
corresponding component must be antisymmetric with respect to q-Pi- 

If the system is parametrized in its rest frame according to Eqs. (A. 4)- 
(A.6), the scalar product q-Pi is given by 

q-Pi=qz q ^p] = iqz q M i , (D.6) 

and the (anti)symmetry of the components in q-Pi is mapped onto an 
(anti) symmetry with respect to the angular variable z q . By additionally 
choosing the momentum partitioning rj± = 0.5 (cf. Sec. A. 2), the compo- 
nents become either real or purely imaginary (see, e.g. [18]). 

These properties allow to restrict the homogeneous BSE in such a way 
that only states of one C-parity appear as solutions. In most cases, this is 
done by expanding the components into Chebyshev polynomials of even or 
odd order depending on the desired symmetry properties (see, e.g., [28, 40, 
51]). 

In the case of the vertex BSE, the inhomogeneous BSA has the same 
symmetry properties as the homogeneous, with the addition that only states 
which have overlap with the inhomogeneous term Tq can contribute, as 
exemplified in Sec. 5.4.2. Thus, if the momentum partitioning r/± = 0.5, 
only poles corresponding to homogeneous amplitudes which have the same 
z-symmetry as Tq appear in the inhomogeneous BSA. 
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Dependence of meson masses 
on cj 

In Figs. E.l - E.6, we show the dependence of all states investigated in 
Sec. 5.3 on the MT parameter u, defined via Eq. (B.2), in comparison to 
the experimental values of the masses. The experimental data is repre- 
sented by black dotted lines, the corresponding errors are indicated by the 
grey bands. The symbols without errorbars are obtained by direct calcula- 
tions, i.e. the masses are within the numerically accessible region, while the 
symbols with errorbars are the results of extrapolations of the eigenvalue 
curves, as described in Sec. 5.2.3. 

The solid lines are linear fits in u to the data points, weighted by the 
appropriate error. In the case of the direct caclulation, the error is deter- 
mined by the numerical precision of the eigenvalue, and estimated by 10~ 7 . 
We consider a total of three states (ground state and two excitations) for 
each J PC for all quark masses listed in App. B, except for J PC = 1 where 
three excitations are calculated. 

In the isoscalar case, six states are shown (eight for J PC = 1 ), which 
correspond to the pure ss states (symbols ( v ) , □, A, <) and the pure nil states 
(symbols o, V, o, >), since we do not consider mixing. 
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Figure E.l: J PC = + , see text. 



84 



APPENDIX E. DEPENDENCE OF MESON MASSES ON u 



Isovector 



Isoscalar 



M[GeV] 




p(770)° 



M[GeV] 



p(1700)° 2.0 



p(1450) u 



1.5 



1.0 



0.3 0.4 0.5 0.6 



w[GeW] 0.5 




w(2205) 
0(2175) 



w(1960) 



0(1680) 
tot 1650) 



£i>(1420) 



0(1020) 



O o w(782) 



0.3 0.4 0.5 0.6 



tj[GeV] 



Charmonium 




Bottomonium 



0(3770) 
0(2S) 



WIS) 



M[GeV] 



10.2- 



10.0- 



0.3 0.4 0.5 0.6 



to[GeV] 




Y(3S) 



Y(2S) 



Y(1S) 



0.3 0.4 0.5 0.6 



to[GeV] 



Figure E.2: J PC = 1 , see text. 
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Figure E.3: J PC = ++ , see text. In the isoscalar channel, the experimental 
error of /o(600) is shaded in lighter gray, since it overlaps with the error of 
/ (980). 
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Figure E.4: J PC = 1 ++ , see text. 
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Isovector Isoscalar 
M[GeV] M[GeV] 




Figure E.6: J PC = 2++, see text. The mass of 72(1430) (isoscalar) is given 
in [86] as M ~ 1430 MeV, and no experimental error is reported. This is 
indicated by omitting the dotted line. 
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